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1.  INTRODUCTION 


Many  years  ago  a  computer  program  was  developed  at  NADC  (Huang  1969)  which 
computes  the  steady  state  shape  of  the  TACAMO  towline  assuming  an  aircraft  orbiting  in  a 
circular  orbit  at  constant  altitude  and  speed.  Although  the  towline  configuration  was  never 
checked  experimentally,  the  vertical  separation  between  towplane  and  conical  drogue  at  the 
bottom  of  the  towline  did  compare  favorably  with  a  limited  number  of  measurements.  Bickel 
et  al  (1971)  using  this  code  along  with  a  VLF  propagation  program  (Pappert  and  Shockey 
1971)  presented  some  calculations  for  the  vertical  electric  field  at  the  ground  produced  by  a 
TACAMO  antenna  and  compared  those  results  with  fields  generated  by  an  elevated  rotating 
dipole.  To  the  authors’  knowledge  those  are  the  only  long  range  propagation  comparisons 
made  with  a  rationally  determined  TACAMO  antenna. 

The  VLF  propagation  code  mentioned  above  was  developed  for  calculating  VLF 
fields  produced  by  antennas  of  arbitrary  length,  shape  and  elevation.  The  program  was  based 
on  simple  segmentation  of  the  antenna  with  each  segment  acting  as  a  short  dipole.  The  field 
strength  at  any  distance  from  the  antenna  was  then  calculated  as  the  phasor  sum  of  the  contri¬ 
butions  from  each  segment.  To  perform  the  calculations  presented  by  Bickel  et  al  (1971),  the 
positions  and  orientations  of  the  antenna  segments  had  to  be  manually  input  to  the  propaga¬ 
tion  code  along  with  the  current  moment  for  each  segment.  The  latter  required  an  assumed 
current  amplitude  and  current  distribution  along  the  antenna.  In  particular,  no  rationale  was 
used  to  relate  current  amplitude  to  radiated  power.  The  latter  is  generally  the  quantity  assumed 
in  systems  calculations. 

The  purpose  of  the  present  report  is  to  describe  a  computer  program  called  TWIRE 
which  combines  several  programs  and  facilitates  calculations  of  the  type  made  by  Bickel  et  al 
(1971)  and  which  uses  as  one  of  its  fundamental  inputs  the  power  radiated  by  the  antenna. 
Specifically,  TWIRE  combines  the  configuration  code  of  NADC  with  a  simple  radiation  resis¬ 
tance  program  (Pappert  1986)  and  a  fast  mode  conversion  propagation  code  (Ferguson  & 
Snyder,  1980).  TWIRE  also  makes  allowance  for  the  counterpoise  which  is  typically  about  an 
eighteenth  of  the  rf  wavelength  (Fern  1986). 

As  mentioned  above,  the  NADC  code  calculates  the  steady  state  antenna  configuration 
by  assuming  an  aircraft  executing  a  circular  orbit  at  constant  altitude  and  speed.  The  radiation 
resistance  code  assumes  a  sinusoidal  current  distribution  along  the  antenna  and  is  strictly  valid 
for  a  flat  infinitely  conducting  ground.  Nevertheless,  it  is  believed  to  be  a  reasonable  approxi¬ 
mation  for  high  conductivity  ground  (e.g.,  seawater).  Towline  plus  counterpoise  is  assumed  to 
be  one-half  of  the  rf  wavelength.  Propagation  in  either  laterally  homogeneous  or  inhomogene¬ 
ous  environments  is  allowed  for  by  the  fast-mode  conversion  program. 

The  NADC  program  determines  steady  state  configurations  of  the  towline  from  force 
equations  by  an  iteration  procedure  which  uses  as  starting  conditions  the  radius  and  altitude  of 
the  conical  drogue  at  the  bottom  of  the  towline.  To  achieve  convergence  it  is  generally  neces¬ 
sary  to  begin  with  good  drogue  starting  conditions.  That  requirement,  unfortunately,  prevents 
full  automation  of  TWIRE.  An  effective  operating  procedure  is  to  determine  starting  condi¬ 
tions  by  running  the  NADC  code  separately  for  a  set  of  operating  conditions,  and  to  then 
employ  those  starting  conditions  in  TWIRE. 

A  brief  description  of  the  NADC  model  is  given  in  the  following  section.  The  radia¬ 
tion  resistance  formula  is  summarized  in  section  3  and  formulas  for  mode  sums  are  given  in 
section  4.  A  brief  program  description  is  presented  in  section  5  and  sample  input  and  output  is 
discussed  in  section  6.  Sample  results  are  given  in  section  7  and  a  summary  including  areas  of 
improvement  concludes  the  study  in  section  8.  A  listing  of  TWIRE  is  given  in  appendix  A  and 
for  the  purpose  of  running  separately  to  obtain  good  drogue  starting  conditions,  the  NADC 
program  is  listed  in  appendix  B. 
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2.  BRIEF  DESCRIPTION  OF  THE  NADC  MODEL 


Several  of  the  assumptions  used  in  constructing  the  NADC  model  are  (Huang  1969): 

1.  Towline  is  inextensible  but  perfectly  flexible. 

2.  Towplane  travels  in  a  perfectly  circular  path  at  a  constant  altitude  and  speed  with 
no  towline  pay-out  or  reel-in. 

3.  No  wind  or  no  wind  shear  present. 

4.  Air  density  varies  with  altitude  as  in  a  standard  atmosphere. 

Working  with  a  cylindrical  coordinate  system  Huang  (1969)  derived  the  steady  state 
configuration  of  the  towline  by  summing  up  the  forces  acting  on  an  element  of  the  towline  AS 
and  equating  the  sum  to  zero.  His  vector  force  equation  is; 

at  d  , .  , .  d  ,  -  , ,  _  02? 

CdP  -  |V_l|V^  +  rrCfP  -  IVReilVR^,  -  ^gk  -  p  -^  =  0  ,  (I) 

where 

s  =  distance  along  towline  measured  from  its  bottom. 

T  =  towline  tension  at  s. 

d  =  diameter  of  towline, 

p  =  air  density. 

Vj_  =  component  of  the  relative  air  velocity  perpendicular  to  the  towline. 

^  Rel  =  relative  air  velocity. 

Cjj  -  drag  coefficient  (dimensionless). 

Cf  =  skin  friction  coefficient  (dimensionless), 

p  =  towline  mass/ unit  length, 

g  =  gravitational  acceleration, 

ic  =  unit  vector  in  z  (altitude)  direction. 

r  =  r  u^  =  radius  vector  from  axis  of  rotation  to  point  on  towline, 

u^  =  unit  vector  from  axis  of  rotation  to  point  on  towline, 

r  =  towline  radius  measured  from  axis  of  rotation. 

When  reduced  to  component  form  equation  (I)  becomes  (Huang  1969) 
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(2) 


(TRT  ^  Tr(0')'  +  Cd  y  p(l  -  (reo-)'  ~T^r’w~Q’  +  ^  o  . 

(TreO'  +  Tr'0'  rrCf  y  p(rto)2  -  Cp  ^  p(ra>)2(I  -  (r0')2)-''  -  =  0  .  (3) 

(TzT  +  Cd  y  pr-^aj20'z'(l  -  (r0')2)'  2  _  =  0  .  (4) 

A  derivative  with  respect  to  s  is  indicated  by  the  prime.  0  is  the  polar  angle,  z  the 

altitude  of  the  towline  at  position  s  and  w  is  the  angular  speed  of  the  orbiting  towplane.  Huang 

(1969)  solved  these  equations  by  an  iteration  procedure  beginning  at  the  drogue  end.  As  men¬ 
tioned  in  the  introduction,  this  generally  requires  good  drogue  starting  conditions.  To  complete 
the  description  of  the  problem,  force  equations  and  parameters  appropriate  to  the  drogue  are 
also  required.  Unfortunately,  the  NADC  report  does  not  make  it  clear  how  the  program  could 
be  altered  to  accommodate  other  drogues,  nor  is  there  sufficient  discussion  to  know  where 
many  of  the  drogue  parameters  came  from  (presumably  some  were  determined  experimen¬ 
tally).  Whether  or  not  the  equations  and  parameters  used  in  the  NADC  report  would  be  suit¬ 
able,  for  example,  to  airborne  command  post  towlines  and  drogues,  is  not  known. 

Under  some  operating  conditions,  it  is  known  that  two  steady  state  solutions  to 
equations  (2)  through  (4)  exist  (Huang  1969  and  Fern  1986).  For  the  conditions  used  as  exam¬ 
ples  in  the  results  section  of  this  report,  it  is  believed  that  only  a  single  solution  exists.  An  often 
observed  operational  behavior  is  a  yo-yo  effect  associated  with  an  orbiting  aircraft.  Fern  (1986) 
gives  an  insightful  discussion  of  the  effect’s  origin,  pointing  out  how  the  lift  on  the  towline 
varies  from  upwind  to  downwind  portions  of  the  aircraft  orbit  and  how  quite  large  excursions 
of  the  towline  could  occur  if  the  operating  conditions  are  such  that  two  steady  state  configura¬ 
tions  exist  simultaneously.  Fern  (1986),  using  the  electromagnetic  code  NEC  (Burke  and 
Poggio  1977),  also  studied  electrical  properties  of  the  TACAMO  antenna.  Although  superior 
to  the  radiation  resistance  calculation  used  in  this  report,  NEC  would  probably  be  too  cumber¬ 
some  to  use  in  combination  with  the  fast  mode  conversion  program.  Perhaps,  though, 
MININEC  (Rockway  and  Logan  1986)  could  be  incorporated  in  the  future. 


3.  RADIATION  RESISTANCE 


This  section  contains  the  formula  for  the  radiation  resistance  used  in  this  report  for 
thin  antennas  of  arbitrary  elevation  and  orientation  over  perfectly  conducting  ground.  The 
formula  is  based  on  segmentation  of  the  antenna  and  figure  1  shows  the  configuration  for  the 
n'^  dipole  of  current  moment  The  dipole  is  located  at  (x^,  y^,  z^)  with  orientation  <f>^  and 
7n  relative  to  the  x  and  z  axes,  respectively.  For  the  special  case  of  a  perfectly  conducting 
ground  the  time  averaged  radiated  power  is  (see  reference  8): 
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Xn 

-  Xr„  ■ 

\ 
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Also,  k  is  the  free  space  wavenumber  and  the  J’s  are  half  integer  order  Bessel  functions 
of  the  first  kind  which  are  expressible  in  terms  of  sines  and  cosines  and,  therefore,  easily  calcu¬ 
lated.  In  the  case  of  a  single  dipole,  equation  (5)  is  reduced  to  the  result  (the  dipole  subscript  is 
omitted  in  (8)  and  (9)  so  that  7  is  the  dipole  orientation  relative  to  the  z  axis): 


Pw  =  20  k2  m2  f(z\7)  . 

where  z*  =  2kz  with  z  the  dipole  height  above  ground  and  where 
f(z*.7)  =  I  + - ^(sin(z*)  -  z*cosz*)  COS27 

L  J 

+  11+  — r  (( I  -  (z*)2)  sinz*  z*cosz*)  si 
|_  2(z*)2  'j 


(8) 


sin27 


(9) 


To  complete  the  description  between  radiated  power  and  antenna  current,  the  current 
distribution  is  assumed  sinusoidal,  a  reasonable  assumption  for  thin  antennas.  In  particular,  as 
mentioned  in  the  introduction,  the  towline  length  plus  counterpoise  is  taken  to  be  one-half  of 
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the  free  space  rf  wavelength.  That  distribution  coupled  w'ith  equation  (5)  then  suffices  to 
determine  the  current  amplitude  in  terms  of  the  radiated  pow'er. 


4.  MODE  SUMS 

In  the  following  (x,  z)  is  the  plane  of  propagation  with  x  the  range  coordinate  and  posi¬ 
tive  z  directed  towards  the  ionosphere  with  z  =  0  the  ground.  The  x  coordinate  of  the  towplane 
is  taken  to  be  x  =  0.  For  a  slab  mode  conversion  model  with  x„  the  beginning  of  the  n'^  slab 
and  with  the  first  slab  described  by  the  region  x  <  X2  the  mode  sums  for  a  W  segmented 
antenna  may  be  written  as 


W 


E'n"  =  QXC(zR)exp(-ikSn'x)  X 


(Zw)s'n7„  +  A.3<„'^f3<h  (z^)sin7^cosd>^  j  , 


(10) 


E'P»  =  Q  X  a 


(p) 


J  m 

w 


;(p) 


f^P>  (ZR)exp(-ik(Sf  )X2  +  S]P>  (X  -  Xp)) 


Mu,exp(ikS*Ox^)  f 

X  -r - ^ - — -.-n-"  (^w)  cosy^  +  fjO*  (z„)  siny^sind)^. 


w=I 


(x  -  X 
sin  — 


J\‘ 


J  sin7^C0Sd.„J  ,  p^\  , 


(11) 


where  j  and  m  are  mode  indices  and  w  is  the  segment  index.  The  barred  quantities  represent 
midpoint  values  for  the  w'*’  segment.  The  superscripts  denote  slab  number  and  the  subscript,  n. 
denotes  the  electric  field  component  at  the  receiver  with  the  convention,  n  =  1  —  EZ,  n  =  2  -* 
EY,  n  =  3  ^  EX.  The  A.’s  and  fs  are  excitation  factors  and  height  gains,  respectively.  The  first 
subscripted  index  on  X  pertains  to  the  orientation  of  the  transmitter  with  X,  being  the  vertical 
dipole  excitation  factor,  Xj  the  broadside  excitation  factor  and  X3  the  end-on  excitation 
factor.  Similarly,  f,  is  the  modal  height  gain  for  EZ,  f2  the  modal  height  gain  for  EY  and  f3 
the  modal  height  gain  for  EX.  With  allowance  for  a  different  indexing  convention  and  height 
gain  normalization,  formulas  for  the  X’s  and  the  fs  will  be  found  in  reference  9.  The  aj^’  are 
cumulative  mode  conversion  coefficients  and  physically  represent  the  accumulative  conversion 
from  a  unit  amplitude  wave  in  mode  m  in  the  transmitter  region  to  mode  j  in  the  p'*’  slab. 
Reference  1 1  shows  how  the  cumulative  conversion  coefficients  are  calculated  by  the  “Fast 
MC”  method.  Other  quantities  appearing  in  equations  (10)  and  (I  I)  are  the  dipole  moment.  M; 
the  earth’s  radius,  a;  the  sine  of  the  modal  eigenangle,  S;  the  free  space  wavenumber,  k;  the 
dipole  segment  orientation  factors,  (7,  d>},  defined  in  the  previous  section;  the  receiver  altitude, 
Zr  and  i  =  \/~T.  The  excitation  factors  and  eigenangle  inputs  are  obtained  from  the  computer 
code  MODESRCH  developed  by  Morfitt  and  Shellman  (1976).  With  the  dipole  moment 
expressed  in  units  of  amp-meters  and  the  electric  field  strength  E  in  microvolts/  m,  the  constant 
Q  is 


Q  =  2.853  X  10  2 

where  f  is  the  frequency  in  kilohertz. 


(12) 


Equation  (10)  is  used  for  laterally  homogeneous  waveguide  caicviations  and  equation 
(11)  for  laterally  inhomogeneous  waveguide  calculations.  Sample  results  are  given  in  a  later 
section.  An  additional  point  is  that  off  axis  antenna  effects  have  been  ignored  in  equations  (10) 
and  (11).  That  should  be  a  good  approximation  for  ranges  of  several  hundred  kilometers  or 
greater. 


5.  BRIEF  PROGRAM  DESCRIPTION 

Appendix  A  contains  a  listing  of  TWIRE.  The  nucleus  of  this  code  is  the  FASTMC 
code  (Ferguson  and  Snyder  1980)  and  in  this  section  only  those  program  elements  of  FASTMC 
which  have  been  altered  or  program  elements  which  have  been  added  to  FASTMC  will  be 
discussed. 

MAIN  controls  the  program  flow.  Data  is  supplied  to  MAIN  via  namelist/ datum  . 

The  latter  is  identical  to  its  counterpart  in  FASTMC.  However,  the  quantities  al.  inch  theta 
and  talt  which  appear  in  namelist  are  not  used  in  TWIRE.  Main  calls  TACAMO,  COORD 
and  RPOWER  which  are  not  program  elements  of  FASTMC. 

TACAMO  is  the  NADC  program  for  calculating  the  steady  state  towline  configuration 
for  an  aircraft  orbiting  in  a  circle  at  constant  altitude  and  speed.  Data  is  supplied  to  TACAMO 
via  namelist  tacin  .  The  namelist  inputs  are; 

v  =  aircraft  speed  in  knots. 

rpl  =  aircraft  radius  in  ft. 

zpl  =  aircraft  altitude  in  ft. 

rzd  =  drogue  starting  radius  in  ft.  If  not  sufficiently  close  to  final  iterate  the 
program  will  abort. 

zzd  =  drogue  starting  altitude  in  ft.  If  not  sufficiently  close  to  final  iterate  the 
program  will  abort. 

amg  =  drogue  weight  in  lbs  (typically  100  lbs). 

psi  =  rotation  angle  (degrees)  discussed  subsequently  in  connection  with  the 
subroutine  coord. 

iclock  =  0  if  aircraft  executes  a  counterclockwise  orbit,  #  0  otherwise. 

itrset  =  number  of  times  Huang’s  original  iteration  scheme  is  used  before  changing 
to  a  two  dimensional  Newton  iteration  which  has  been  added  to  the  original 
NADC  program.  Generally,  itrset  has  been  set  to  zero  in  carrying  out  the 
sample  calculations  in  this  study. 
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cpl  =  counterpoise  length  in  ft. 

chicpl  =  inclination  of  counterpoise  in  degrees  from  horizontal  (90°  <  chicpl  <  0° 
with  0°  horizontal  and  90°  vertical). 

Aerodynamic  characteristics  of  the  conical  drogue,  determined  by  parameters  not 
listed  above,  are  fixed  in  data  statements  and  were  apparently  derived  from  data  presented  in 
reference  7.  Towline  drag  and  mass  coefficients  are  also  fixed  in  data  statements.  Lengths  of 
the  towline  segments  in  ft.  are  given  in  the  data  statement  dstore  for  segments  1  through  1 7  and 
are  given  by  dds  for  segments  greater  than  17.  These  values  can  be  altered  to  accommodate 
individual  needs. 

The  Newton  iteration  mentioned  above  is  developed  from  the  definitions 


f|  =  Z|(Zd,rd)  -  Zpi  , 

(13) 

f,  =  r,(Zd,rd)  -  rp,  , 

(14) 

where  z,  is  the  iterated  plane  altitude  and  r,  the  iterated  plane  radius  corresponding  to  the 
starting  drogue  conditions  z^  and  r^.  Also,  Zpj  and  fp]  are  the  true  plane  altitude  and  radius 
respectively.  Expanding  f,  and  {2  about  Zj  and  r^  yields 


9zi  9Z| 

+  gzj  +  —  Srj  «  0  , 

(15) 

9r|  9r| 

9Zd  orj 

(16) 

Solution  of  equations  (15)  and  (16)  gives 


9r 


9z 


A  = 


■d  '^'d 

3Z(  9r,  9r,  9z, 


I 


9zd  9rd  9zd  9rd 


(17) 

(18) 


TACAMO  is  also  included  as  a  separate  entry  in  Appendix  B.  As  mentioned  previ¬ 
ously,  it  is  best  when  working  with  new  operating  conditions  to  run  TACAMO  by  itself  to  get 
starting  conditions  for  TWIRE.  A  possible  step  to  fuller  automation  would  be  to  develop  a 
multidimensional  grid  of  starting  values  for  a  cross  section  of  operating  conditions  and  to  use 
interpolation  to  determine  drogue  starting  values  for  intermediate  operating  conditions.  In  that 
connection,  operating  frequencies  range  from  about  17  kHz  to  30  kHz  and  TACAMO  aircraft 
generally  operate  at  air  speeds  between  about  170  and  250  knots  at  altitudes  between  about 
20,000  and  40,000  ft.  with  radii  >  4,000  ft.  (Fern  1986). 

COORD  performs  the  following  coordinate  transformations.  In  the  schematic  below, 
the  towplane's  (x,y)  coordinates  are  denoted  by  (Xp,yp).  The  axis  are  first  rotated  by  the  angle 
6  =  tan  '  (yp/Xp),  so  that  in  the  x',  y'  system  y'  =  0. 

The  y'  axis  is  the  direction  of  the  towplane  for  counterclockwise  rotation  (i.e., 
iclock  =  0).  Coordinates  of  the  point  (x,y)  in  the  unprimed  system  become  in  the  prime  system 
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XCOS0  +  ysind  , 

(19) 

-xsinfl  +  ycosS  . 

(20) 

The  prime  coordinate  system  is  next  translated  so  that  the  towplane  is  at  the  origin.  That  is 

XT  =  x'-Xp,  (21) 

y-f  =  y'  ■  (22) 

If  the  airplane  is  moving  in  a  clockwise  rotation  (iclock  =  1),  the  y'  coordinates  are 
reversed  in  sign.  Let  ijt  (psi)  be  the  direction  between  the  x'  direction  (i.e.,  direction  from  center 
of  orbit  to  towplane)  and  the  direction  of  propagation,  denoted  by  x". 

Then, 


x"  =  xjcosilt  +  y  jsini/> ,  (23) 

y"  =  -x-f  cosi/r  +  y  jcosi/r .  (24) 

Values  of  ij/  at  cardinal  points  on  the  orbit  with  x"  directed  to  the  right  are  indicated 

below. 
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The  counterpoise  coordinates  are  also  calculated  in  COORD  subject  to  the  assumption 
that  the  counterpoise  is  in  the  direction  -v/  |vj,  where  v  is  the  towplane’s  velocity  and  droops  at 
an  angle  chicp  with  the  horizontal. 

RPOWER  is  the  program  element  where  equation  (5)  is  implemented  subject  to  the 
assumption  that  the  antenna  length  plus  counterpoise  is  one-half  the  rf  wavelength  and  that  the 
current  is  sinusoidally  distributed.  Outputs  of  the  subroutine  are  the  current  moments  (idl)  of 
the  antenna  segments  for  a  given  radiated  power  input,  pin. 

MCFLDS  is  an  original  program  element  of  the  program  FASTMC  which  in 
TWIRE  has  been  modified  to  allow  for  antenna  segmentation.  IF  iopt  =  1  (range  option  flag), 
MCFLDS  is  called  from  MAIN  for  each  range  between  dmin  and  dmax  at  deltad  intervals 
and  the  electric  field  components  EZ  or  EY  are  computed  according  to  equations  (10)  or  (1 1). 
Field  values  in  complex  form  are  contained  in  SUMOUT  and  field  amplitudes  in  Amp.  If  iopt 
=  2  (height  gain  option),  MCFLDS  is  used  as  an  intermediate  step  in  calculating  the  electric 
field  components.  Final  field  calculations  in  this  instance  are  performed  in  MAIN.  Again,  field 
values  in  complex  form  are  contained  in  SUMOUT  and  field  amplitudes  in  Amp. 


6.  SAMPLE  INPUT  AND  OUTPUT 

Three  files  are  required  for  input  to  TWIRE.  The  three  files  are  called  twzy.in, 
twzy.com  and  cards. 30.  They  represent  input  for  a  range  calculation  of  the  vertical  electric 
field  component  EZ  and  the  horizontal  electric  field  component  EY  for  a  laterally  homogene¬ 
ous  guide.  The  file  twzy.in  first  calls  for  the  file  twzy.com  and  then  for  data  contained  in  name- 
list  tacin,  which  is  used  for  the  towline  configuration  calculation.  Listed  in  order  in  namelist 
tacin  are: 

V  =  towplane  velocity  in  knots, 
rpl  =  towplane  radius  in  feet, 
zpl  =  towplane  altitude  in  feet, 
rzd  =  drogue  starting  radius  in  feet, 
zzd  =  drogue  starting  altitude  in  feet, 
amg  =  drogue  weight  in  pounds. 

psi  =  angle  in  degrees  between  radius  vector  from  center  of  towplane’s  orbit  to 
towplane  and  the  direction  of  rf  propagation. 

iclock  =  flag  for  clockwise  (1)  or  counterclockwise  (0)  towplane  rotation. 

chicp  =  droop  angle  of  the  counterpoise  measured  in  degrees  from  the  horizontal. 

cpl  =  counterpoise  length  in  feet. 
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The  namelist  tacin  is  then  repeated  a  second  time  as  the  configuration  calculation  is 
repeated  for  each  component  calculated. 

The  file  twzy.com  begins  with  the  character  string  “name”.  Its  remaining  contents  in 
this  instance  are; 

rcomp  =  elearic  field  component  calculated  (1  for  EZ,  2  for  EY). 
nprint  =  printout  flag 
=  0  none 

=  1  slab  identification  data  and  field  amplitude  and  phase 
=  3,4  intermediate  mode  conversion  data  used  for  diagnostics  only, 
naplot  =  0  amplitude  plot  flag  (0  for  no  plot,  1  for  plot), 
dmin  =  minimum  range  in  Mm. 
dmax  =  maximum  range  in  Mm. 

nrptsl  =  number  of  ranges  calculated  at  equal  increments  between  dmin  and  dmax. 
ampmax  =  upper  limit  of  field  amplitude  plot  dB//iV/m. 
ampmin  =  lower  limit  of  field  amplitude  plot  in  dB/)uV/m. 
ralt  =  receiver  altitude  in  km. 

cards.30  =  waveguide  mode  constant  data  (see  e.g.,  reference  6). 

The  calculations  are  repeated  for  rcomp  =  2  (i.e.,  for  EY  calculations). 

The  waveguide  input  file,  cards.30,  obtained  from  reference  6  begins  with  DATA 
followed  by  a  descriptor  title.  The  next  line  identifies  the  slab  properties  as  follows; 

R  =  range  location  of  slab  in  km. 

F  =  frequency  in  kHz. 

A  =  propagation  azimuth  from  magnetic  north  (degrees). 

C  =  geomagnetic  field  codip  (degrees). 

M  =  geomagnetic  field  strength  (W/m^). 

S  =  ground  conductivity  (S/  m). 

E  =  ground  dielectric  constant. 

T  =  Top  height  used  for  waveguide  calculations  (km). 
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INPUT  FILE: 


Cwzy.com 
&cacin 
v-200. . 
rpl-4000. , 
zpl-20000 . , 
rzd-700.0, 
zzd-9000. , 
amg-100 . , 
psl-180.0, 
lclock>-0 , 
chicp-0 . , 
cpl-1800. , 
&end 

y 

y 

&Cacin 
v-200 . , 
rpl-4000 . , 
zpl-20000. . 
rzd-700 . 0 , 
zzd-9000 . , 
amg-lOO . , 
pst-180.0, 
lclock-0 , 
chlcp-0 . , 
cpl-1800. , 
&end 


TWZY.IN 
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INPUT  FILE:  TWZY.COM 


name 

&datum  rcomp— 1 , nprint— 1 .naplot-O , 
dmin-0 . 000 , dmax-5 . 000 , 
nrptsl-ll, 

aiDpmax~80 . ,  ampmin—O . , 
ralt-9.144,  &en<l 
input  cards . 30 
name 

&daCum  rcomp— 2 , &end 
start 
quit 
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INPUT  FILE;  CARDS. 30 


DATA 

towline  test 

R  .000  F  30.0000  A  90.000  C  40.000  M  .500E-04  S  4.640E+00  E  81.0  T  87.0 

1  89.97300  -6.722872  1 . 14900804E-05-9 . 69784160E-05-5 . 66454815E-13  8 . 91278878E-16 

2  89.97300  -6.722872  4 . 94460073E-09  5 . 55613688E-09  1 . 00000072E+00-7 . 01791265E-08 

1  89.79732  -6.138471  1 . 39372627E-04-4. 28004772E-04-9 . 34471576E- 13- 3 .41234067E- 13 

2  89.79732  -6 . 138471-9 . 26029919E-09-1 . 90273877E-08  1 . 00000119E+00  9 . 79737678E-08 

1  89.91204  -3.159701-5.67914452E-04-1.12545611E-02-4.18753747E-12-1.40717722E-12 

2  89.91204  -3.159701  1 . 34759489E-07  1 . 77823267E-C7  1 . 00000048E+00-4 . 46208873E-07 

1  88.82600  -1.083422  2 . 33682222E-03-5 .44904545E-03-2 . 69803867E- 11-6 .00914040E- 14 

2  88.82600  -1 . 083422-2 . 19763109E-07-3 . 34162451E-07  1 . 00000072E+00  7 . 85546206E-09 

1  84.71127  -. 222291-2. 45107268E-03-1.99295953E-02-1.35681544E-11-9.35714700E-12 

2  84.71127  -.222291  2 . 98093312E-07  4.92028562E-07  1 . 00000048E+00  2 . 00596716E-07 

1  82.74648  -.273032  2 . 91138096E-03-4. 90130810E-03-1. 01852499E-10  9 . 41601918E- 12 

2  82.74648  - .273032-4. 07881686E-07-6.45559965E-07  1.00000012E+00-5.70226248E-08 

1  80.40218  -. 234551-3. 23274592E-03-1.51210632E-02-5.33708217E-11-3.51053041E-11 

2  80.40218  -.234551  5 . 61018396E-07  8 . 20390710E-07  1 . 00000036E+00  1 . 91906668E-08 

1  78.70553  -.293832  3 . 81615874E-03-6 .45386567E-03-1 . 83607657E-10  3 . 90313511E- 11 

2  78.70553  -. 293832-6 . 92323852E-07-9 . 63369075E-07  1 .00000048E+00  2 .79036698E-07 

1  76.58138  - .261791-4. 16685035E-03-1.15461908E-02-1.42259260E-10-8.68795105E-11 

2  76.58138  -.261791  9 . 04359183E-07  1. 10827352E-06  9 . 99999583E-01  3 . 27973055E-07 

1  74.96912  -.345282  4 . 87595052E-03-8 . 77621677E-03-2 . 45324872E- 10  9 . 70813221E-11 

2  74.96912  -. 345282-1. 07347228E-06-1.22331073E-06  9 . 99999821E-01  4.46039110E-07 

1  72.92796  -. 296672-5 . 04938047E-03-8 . 35422613E-03-3 . 03403247E-10-1 . 68299430E-10 

2  72.92796  -.296672  1 . 32530784E-06  1 . 27687861E-06  1 . 00000048E+00  9 . 59014642E-07 

1  71.32272  -.409471  5 . 76619152E-03-1 . 14735551E-02-2 . 64180205E-10  1.84621096E-10 

2  71.32272  - .409471-1. 53844314E-06-1.33109688E-06  1 . 00000048 E-^00  9 . 05900492E-07 

1  69.33108  -. 339462-5. 64039499E-03-5.37072308E-03-5.54849555E-10-2.73569306E-10 

2  69.33108  -.339462  1 . 79977269E-06  1. 25654230E-06  9 . 99999106E-01  1 . 31630159E-07 

1  67.68570  -.479001  6 . 25931984E-03-1.42788989E-02-2 . 23180446E-10  2 . 91375785E- 10 

2  67.68570  - .479001-2. 05986134E-06-1.21624373E-06  9 . 99997795E-01  2 . 21355563E-07 

1  65.73211  - .392672-5. 82792237E-03-2.63220887E-03-9.07945108E-10-3.89095423E-10 

2  65.73211  -.392672  2 . 30049409E-06  1 . 01230273E-06  9 . 99995053E-01-3 .48421486E-06 

1  64.01385  -.546571  6 .22112770E-03-1.70211233E-02-1.11896700E-10  3 . 96509547E- 10 

2  64.01385  -. 546571-2. 59987064E-06-8.40697567E-07  1.00000191E-H00-1.35956998E-05 

1  62.08826  - .459332-5. 56169311E-03-1.49533094E-04-1.37232026E-09-4.94906449E-10 

2  62.08826  -.459332  2 . 79967230E-06  5 . 28230146E-07  1 . 00002074E+00  2 .41414S64E-06 

r  40.0 
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ramax  =  maximum  receiver  altitude  in  km  (used  with  iopt  =  2). 

ratic  =  tic  mark  spacing  on  receiver  altitude  scale  (used  with  iopt  =  2). 

dist  =  range  locations  in  Mm  at  which  height  gains  are  calculated  (used  with 
iopt  =  2). 

nrd  =  numbers  of  dist  (used  with  iopt  =  2). 

nrpts2  =  number  of  height  gain  points  calculated  between  ramin  and  ramax  (used  with 
iopt  =  2). 

Following  the  quantities  contained  in  the  datum  namelist,  comes  the  slab  identification 
data  described  previously  and  the  tacin  namelist  data  also  described  previously  (except  for 
itrset  which  controls  the  iterations  as  described  in  section  5  —  default  is  0).  Apart  from  input 
quantities  discussed  in  connection  with  namelist  tacin,  the  following  appear: 

abase  =  base  area  of  drogue  in  sq.  ft. 

cddrog  =  drogue  drag  coefficient. 

amug  =  towline  weight  lbs/ ft. 

d  =  towline  diameter  in  ft. 

al  =  towline  length  in  ft. 

dds  =  a  spacing  increment  along  the  towline  in  ft. 

picf  =  towline  skin  friction  coefficient. 

cd  =  towline  drag  coefficient. 

All  of  the  above,  except  al,  are  set  in  data  statements  in  the  subroutine  TACAMO.  The 
quantity  al  is  calculated  subject  to  the  condition  that  the  towline  length  plus  counterpoise  is 
one-half  a  rf  wavelength. 

Following  the  slab  identification  card  is  the  mode  data.  Each  eigenangle  in  degrees 
appears  twice  along  with  four  complex  quantities  from  which  all  excitation  factors  are  calcu¬ 
lated.  The  cards.30  file  is  terminated  with  the  r  =  40.0  character  string.  If  more  slabs  were  used 
in  the  modeling  (i.e.,  lateral  inhomogeneity  allowed  for)  then  slab  identification  and  data  for 
each  additional  slab  would  follow  in  order  with  a  space  separating  each  new  slab  entry  begin¬ 
ning  with  the  range  information  (R  etc.)  record. 

Output  corresponding  to  the  preceding  input  is  shown  on  pages  16  through  21.  The 
output  first  echoes  the  namelist  datum.  Quantities  not  specified  in  the  datum  entry  in  twzy.com 
appear  in  the  output  along  with  their  default  values.  In  order,  quantities  not  discussed  in 
connection  with  the  input  are  (also  see  reference  3): 

iopt  =  range  or  height  gain  option  flat  (1  for  range,  2  for  height  gain). 
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al,  incl.  theta,  talt,  TOPHT  =  quantities  not  used  as  input  in  TWIRE. 
npplot  =  phase  plot  flag  (0  no  plot,  I  plot), 
npdiff  =  not  used  in  present  program, 
nrcurv  =  number  of  curves  on  a  plot  (maximum  of  4). 
sizex  =  size  of  x  axis  in  inches, 

sizey  =  size  of  y  axis  in  inches, 

amptic  =  tic  mark  spacing  on  amplitude  axis, 
phsmax  =  phase  axis  maximum  in  degrees, 
phsmin  =  phase  axis  minimum  in  degrees, 
phstic  =  tic  mark  spacing  on  phase  axis, 
dtic  =  tic  mark  spacing  in  Mm  on  range  axis. 

tlong,  tlat,  rbear,  totape  =  not  used  in  TWIRE. 
power  =  radiated  power  in  kw. 

ramin  =  minimum  receiver  altitude  in  km  (used  with  iopt  =  2). 

The  towline  configuration  follows  and  is  described  by: 
s  =  distance  in  feet  along  towline  as  measured  from  its  bottom, 
r  =  radius  of  towline  in  feet  at  s. 

7  =  altitude  of  towline  in  feet  at  s. 

th  =  azimuth  of  towline  in  radians  measured  in  counterclockwise  direction  at  s. 

The  remaining  quantities  t,  rp,  zp  and  rthp  are  of  no  interest  in  the  present  application. 
In  the  output  it  will  be  seen  that  the  program  iterated  twice,  there  being  three  sets  of  s,  r,  z  and 
th  values  given  at  the  bottom  and  top  of  the  towline.  Final  values  calculated  from  the  bottom 
to  the  top  follow.  Next  comes  imid  which  is  the  towline  segment  where  the  current  is  a  maxi¬ 
mum.  The  midpoint  coordinates  are  listed  for  each  segment  after  the  coordinate  transforma¬ 
tions  discussed  in  section  5  have  been  made.  Also  listed  is  gamma  and  phi  for  each  segment 
where  gamma  is  the  inclination  in  degrees  of  the  segment  from  vertical  and  phi  is  its  azimuth  in 
degrees.  Following  the  midpoint  information  is  a  printout  of  the  electric  field  component  being 
calculated  and  the  receiver  altitude.  The  principal  output  then  follows  showing  dist  in  Mm, 
field  amplitude  in  dB/pV/m  and  its  phase  in  degrees.  The  output  then  repeats  itself  for  the 
second  field  component,  EY. 
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raANSPOraiATION 

>t  xmid  ymid  zmid  gamma  phi 

3991.607  -749.322  9011.749  9.200  66.535 

3994.340  -744.185  9046.258  9.736  60.267 

4004.651  -730.949  9134.634  11.189  49.348 
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7.  SAMPLE  RESULTS 


Shown  in  figures  2  through  5  are  drogue  altitude  and  drogue  radius  for  towline  lengths 
ranging  between  14,000  and  30,000  ft.  The  figures  apply  to  a  towplane  altitude  of  30,000  ft  and 
speed  of  200  knots.  Figure  2  is  for  a  towplane  radius,  r^,  of  4000  ft.,  figure  3  for  r^  =  6000  ft., 
figure  4  for  r^  =  8000  ft,  and  figure  5  for  r^  =  30,000  ft.  which  approximates  straight  flight. 
These  curves  can  be  used  to  obtain  starting  conditions  for  other  operating  conditions  by  grad¬ 
ually  evolving  from  the  conditions  of  the  figures  to  the  desired  operating  conditions  (the  proce¬ 
dure  can  be  very  tedious  in  some  instances).  As  mentioned  previously,  a  possible  step  to  fuller 
automation  of  the  present  program  would  be  to  develop  a  multidimensional  grid  of  starting 
values  for  a  variety  of  operating  conditions,  of  which  figures  2  through  5  are  a  sample,  and  to 
use  interpolation  to  determine  drogue  starting  conditions  for  intermediate  operating  conditions. 

It  will  be  seen  that  there  is  a  rapid  change  in  drogue  radius  and  altitude  which  occurs 
for  a  towline  length  of  about  16,500  ft.  in  figure  3  and  for  a  towline  length  of  about  21,500  ft. 
in  figure  4.  On  the  basis  of  the  examination  made  in  this  study  it  does  not  appear  that  there  are 
two  solutions  in  the  neighborhood  of  those  points  and  the  physical  origin  of  the  rapid  changes 
is  not  known. 

It  is  known  that  a  simple  exponential  ji  -  0.5  km  ',  h'  =  87  km  (notation  of  Wait  & 
Spies,  1964)  profile  does  an  excellent  job  of  predicting  nocturnal  VLF  propagation  to  the  east 
(Papperi  and  Hitney  1988).  Figures  6  through  13  show  results  for  the  vertical  electric  field  at 
the  ground  generated  by  a  T  AC  AMO  antenna  under  several  operating  conditions.  The  results 
are  for  easterly  propagation  (azimuth  =  90°)  and  the  j8  =  0.5  km  ',  h'  =  87  km  profile.  For  all  of 
the  figures  a  seawater  path  is  assumed  (i.e.,  conductivity  =  4.64  s/m,  relative  permittivity  =  81) 
and  the  geomagnetic  field  strength  is  taken  to  be  0.5  Gauss  with  a  dip  of  50°.  Also,  for  all  figures 
the  towplane’s  speed  is  taken  to  be  200  knots.  Figures  6  through  9  are  for  20  kHz  whereas 
figures  10  through  13  are  for  30  kHz.  Results  are  given  for  towplane  altitudes  (Aalt)  of  20,000 
and  30,000  ft.  and  for  towplane  radii  (Rpl)  of  4000  and  8000  ft.  The  counterpoise  length  has 
been  taken  to  be  1800  ft,  at  30  kHz  and  2800  ft.  at  20  kHz.  All  of  the  figures  apply  to  counter¬ 
clockwise  rotation  (Iclock  =  0)  with  i/i  =  0°  (see  section  5)  and  a  radiated  power  (Pwr)  of  1  kW. 
Also  shown  on  the  curves  are  results  for  point  dipole  calculations  with  altitude  and  angular 
orientation  factors  for  the  point  dipole  determined  by  the  towline  segment  which  contains  the 
current  maximum  (i.e.,  imid  discussed  in  section  6).  In  calculating  the  point  dipole  results, 
allowance  has  been  made  for  its  radiation  resistance  dependence  on  height  and  tilt  relative  to 
the  z  axis  as  expressed  by  equation  (9).  Except  in  rare  instances  the  point  dipole  results  agree 
with  the  TACAMO  antenna  calculation  to  better  than  2  dB  and  the  results  are  almost  indistin¬ 
guishable  in  figures  6  and  8.  Interestingly  these  are  the  cases  of  high  verticality  or  low  7’s. 

Although  probably  not  as  much  of  interest  as  the  air-to-ground  transmission  of  EZ  just 
discussed,  figures  14  through  21  show  results  for  the  transverse  electric  (TE)  field  component 
EY  at  30,000  ft  (i.e.,  Ralt  =  30.00).  Other  than  the  receiver  altitude  difference,  the  re.sults  are  for 
the  same  set  of  operating  conditions  which  applied  to  figures  6  through  13.  Shown  again  are 
re.sults  for  the  point  dipole.  It  will  be  seen  that  except  in  rare  instances  the  shapes  of  the  mode 
sum  plots  for  the  TACAMO  antenna  and  point  dipole  calculations  are  very  similar.  Agreement 
in  terms  of  absolute  dB  levels  is  not  quite  as  good  in  this  instance  as  it  was  for  the  EZ  compo¬ 
nent  and  there  are  occurrences  where  the  signal  levels  differ  by  10  dB  or  more. 

Although  the  /3  =  0.5  km  h'  =  87  km  profile  does  not  accurately  predict  westerly 
propagation  (Pappert  and  Hitney  1988),  figures  22  through  25  are  presented  as  illustrating,  for 
that  direction  (azimuth  =  270°)  and  profile,  comparisons  between  the  TACAMO  antenna  and 
point  dipole  calculations.  Figures  22  and  23  show  results  for  Ez  at  the  ground  at  20  kHz  with 
the  towplane  altitude  20,000  ft.  Figure  22  gives  results  for  the  towplane  radius  of  4000  ft.  and 
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figure  23  gives  results  for  the  towplane  radius  of  8000  ft.  The  TACAMO  and  point  dipole 
results  for  the  high  verticality  (low  7)  case  are  in  very  good  agreement.  The  TACAMO  and 
point  dipole  results  for  the  low  verticality  case  (high  7)  shown  in  figure  23  do  not  agree  as  well 
in  this  instance  as  they  did  for  the  corresponding  easterly  propagation  case  (figure  7).  It  is 
known  that  the  VLF  fields  penetrate  to  higher  altitudes  for  westerly  propagation  than  they  do 
for  easterly  propagation.  Thus,  the  modes  are  more  hybrid  for  westerly  propagation  and  that 
may  be  partly  responsible  for  the  larger  differences  noted  for  the  low  verticality  westerly 
propagation. 

Figures  24  and  25  show  results  for  EY  at  30,(X)0  ft.  at  20  kHz  with  the  towplane 
altitude  20,000  ft.  Figure  24  gives  results  for  the  towplane  radius  of  4000  ft.  and  figure  25  gives 
results  for  the  towplane  radius  of  8000  ft.  Again,  the  TACAMO  and  point  dipole  calculations 
for  the  high  verticality  (low  7)  case  are  in  very  good  agreement.  The  differences  between  the 
TACAMO  and  point  dipole  calculations  shown  in  figure  25  for  the  low  verticality  (high  7)  case 
appear  to  be  comparable  to  those  for  the  corresponding  easterly  propagation  case  (figure  15). 

Figures  26  and  27  compare  signal  levels  for  easterly  propagation  for  i/r  =  90°  and  270°. 
The  results  apply  to  a  frequency  of  30  kHz  with  the  towplane  at  20,000  ft.  and  in  an  orbit  of 
radius  4000  ft.  Figure  26  is  for  EZ  at  the  ground  and  figure  27  is  for  EY  at  30,000  ft.  Differ¬ 
ences  between  the  curves  are  measures  of  the  orbital  modulation  expected  as  the  towplane 
executes  its  orbit.  The  modulation  expected  is  one  maximum  and  one  minimum  associated 
with  each  rotation  of  the  towplane.  In  rare  instances,  the  minimum  field  occurs  for  ij/'s  close  to 
0°  and  180°.  Based  on  curves  not  shown  such  a  circumstance  for  the  EY  component  occurs  in 
the  present  example  at  a  range  of  1.425  Mm.  The  expected  modulation  in  this  case  is  twomax- 
ima  and  two  minima  during  the  course  of  one  rotation  of  the  towplane.  Examples  of  each  type 
of  modulation  are  shown  in  figure  28.  The  amplitude  excursions  at  2.6  Mm  for  Ez  and  at  2 
Mm  for  EY  are  consistent  with  expectations  based  on  figures  26  and  27  respectively. 

1  he  preceding  results  have  all  been  for  laterally  uniform  waveguides.  Results  for  a 
laterally  nonuniform  waveguide  are  shown  in  figures  29  and  30.  The  guide  is  characterized  by 
waveguide  modes  corresponding  to  an  azimuth  of  90°  out  to  2  Mm,  between  2  Mm  and  4  Mm 
the  waveguide  modes  are  taken  to  be  those  corresponding  to  an  azimuth  of  270°  and  beyond 
4  Mm  the  waveguide  modes  are  again  taken  to  be  those  characterized  by  an  azimuth  of  90°. 
The  electron  density  in  each  slab  is  described  by  the  P  =  .5km~',  h'  =  87  km  profile  used 
previously.  Also,  the  geomagnetic  field  magnitude  and  dip  as  well  as  the  ground  parameters 
are  identical  to  those  used  previously.  Although  such  a  waveguide  is  physically  unrealistic  it 
can  be  used  as  a  check  on  the  performance  of  the  program  for  a  laterally  inhomogeneous 
environment.  Figure  29  provides  results  at  30  kHz  for  EZ  at  30,000  ft.  The  towplane’s  altitude 
is  30,000  ft.,  its  radius  is  4000  ft.  and  its  speed  is  200  knots.  Figure  30  gives  results  at  30  kHz 
for  the  EY  component  at  30,000  ft.  for  the  same  towplane  conditions.  Also,  shown  on  the  plots 
are  the  point  dipole  results.  The  good  agreement  between  the  TACAMO  antenna  results  and 
the  point  dipole  results  serves  as  a  check  on  the  performance  of  TWIRE  in  a  laterally  inhomo¬ 
geneous  environment. 
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8.  SUMMARY  AND  AREAS  OF  IMPROVEMENT 


A  computer  program  called  TWIRE,  which  combines  a  TACAMO  antenna  configur¬ 
ation  code  (Huang  1969),  a  radiation  resistance  code  (Pappert  1986)  and  a  VLF  propagation 
code  (Ferguson  and  Snyder  1980),  has  been  assembled  for  the  purpose  of  conducting  case 
studies  for  VLF  airborne  transmissions  and  for  the  purpose  of  pointing  out  possible  areas 
of  improvement  for  calculating  such  transmissions.  In  particular,  by  combining  the  three 
codes  TWIRE  provides  a  tool  for  conducting  VLF  airborne  transmission  studies  based  on  a 
rationally  determined  antenna  configuration.  A  major  deficiency  of  the  program  is  that  it  is 
not  fully  automated  because  the  configuration  code  requires  initial  inputs  which  often  have  to 
be  determined  in  advance.  A  possible  scheme  for  overcoming  that  deficiency  would  be  to 
provide  a  multidimensional  grid  of  starting  conditions  for  different  operating  conditions  (i.e., 
towline  length  and  towplane’s  altitude,  speed  and  radius)  and  to  use  an  interpolation  scheme 
for  intermediate  operating  conditions  One  problem  that  would  have  to  be  resolved  in  connec¬ 
tion  with  this  approach  would  be  the  resolution  of  how  to  handle  operating  conditions  when 
more  than  one  stable  towline  configuration  exists.  Other  areas  of  improvement  relate  to  a  fuller 
treatment  of  the  dynamics  of  the  towline.  The  possibility  of  extending  the  steady  state  NADC 
mood  to  include  transient  and  wind  effects  is  such  an  area.  Also,  the  NADC  report  does  not 
fully  describe  the  origin  of  the  drogue  parameters  and  formulas  which  provide  crucial  starting 
conditions  for  the  integrations  necessary  to  determine  the  towline  configuration.  It  is  not  clear, 
therefore,  how  to  adapt  the  NADC  program  to  airborne  systems  other  than  TACAMO.  Better 
documentation  of  the  NADC  code  would,  therefore,  be  of  value. 

Because  of  simplifications  used  in  the  development  of  the  radiation  resistance  calcu¬ 
lations,  TWIRE  is  strictly  valid  for  infinite  ground  conductivity.  This  is  probably  the  easiest 
element  to  improve.  It  is  quite  likely  that  approximate  allowance  for  finite  ground  conductivity 
could  be  made  without  increasing  cpu  time.  It  is  even  possible  that  MININEC  could  be  incor¬ 
porated  in  place  of  the  program  element  RPOWER. 

In  the  present  study,  the  principle  objective  has  been  to  illustrate  the  feasibility  of 
performing  VLF  airborne  transmissions  for  a  rationally  determined  antenna  configuration 
without  specific  regard  to  cpu  time.  It  is  quite  possible  that  TWIRE  can  be  speeded  up  by 
simply  rearranging  do  loops  involving  sums  over  modes  and  antenna  segments. 

Perhaps  the  greatest  utility  of  TWIRE  will  be  to  simply  assess  the  adequacy,  for  field 
calculation  purposes,  of  point  dipoles  as  approximations  to  TACAMO  configurations.  For 
example,  limited  results  given  in  section  7  would  indicate  that  a  sensibly  chosen  point  dipole 
would  suffice  for  many,  if  not  most,  systems  applications. 
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Figure  5.  Towline  length  vs.  drogue  starting 
conditions. 
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Figure  6.  EZ  field  strength  comparisons  between 

TWIRE"  ( - )  and  point  dipole  ( - )  results 

for  an  azimuth  of  90*.  Point  dipde  orientation  and 
altitude  given  by  7  =  17.84“,  <l>  =  -58.69*  and  2  = 
12,945  ft. 
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Figure  8.  EZ  field  strength  comparisons  between 

'TWIRE"  ( - )  and  point  dipole  (•  •  ■  ■)  results 

for  an  azimuth  of  90°.  Point  dipole  orientation  and 
-  altitude  given  by  7  =  19.00°,  (j>  =  -47.18°  and  z  = 
22,627  ft. 
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Figure  11.  EZ  field  strength  comparisons  between 

"TWIRE"  ( - )  and  point  dipole  (•  •■  •)  results 

for  an  azimuth  of  90°.  Point  dipole  orientation  and 
altitude  given  by  7  =  73.08°,  <i>  =  40.75°  and  z  = 
18,287  ft. 


Figure  12.  EZ  field  strength  comparisons  between 

"TWIRE"  ( - )  and  point  dipole  ( - )  results 

for  an  azimuth  of  90°.  Point  dipole  orientation  and 
altitude  given  by  7  =  36.40°,  <l>  =  -21.89°  and  z  = 
26,294  ft. 
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Figure  17.  EY  field  strength  comparisons  between 

'TWIRE"  ( - )  and  point  dipole  (■  •  •  •)  results 

for  an  azimuth  of  90®.  Point  dipole  orientation  and 
altitude  given  by  t  =  55.89°,  <i>  =  -3.63°  and  z  = 
26.037  ft. 
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Figure  18.  EY  field  strength  comparisons  between 

TWIRE"  ( - )  and  point  dipole  (•  •  •  •)  results 

for  an  azimuth  of  90°.  Point  dipole  orientation  and 
altitude  given  by  7  =  38.38°,  ^  =  -29.96°  and  z  = 
16,757  ft. 
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Figure  19.  EY  field  strength  comparisons  between 

"TWIRE"  ( - )  and  point  dipole  (■  •  •  •)  results 

for  an  azimuth  of  90°.  Point  dipole  orientation  and 
altitude  given  by  7  =  73.08°,  (j>  =  40.75°  and  2  = 
18,287  ft. 


Figure  20.  EY  field  strength  comparisons  between 

"TWIRE"  ( - )  and  point  dipole  (•  ■  •  ■)  results 

for  an  azimuth  of  90°.  Point  dipole  orientation  and 
altitude  given  by  7  =  36.40°,  ({>  =  -21.89°  and  z  = 
26,294  ft. 
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Figure  23.  EZ  field  strength  comparisons  between 

"TWIRE"  ( - )  and  point  dipole  (•  •  ■  •)  results 

for  an  azimuth  of  270°.  Point  dipole  orientation 
and  altitude  given  by  ^  =  70.33°,  <i>  =  7.23°  and 
z  =  17,250  ft. 
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Figure  24.  EY  field  strength  comparisons  between 

'TWIRE"  ( - )  and  point  dipole  (•  •  •  •)  results 

for  an  azimuth  of  270°.  Point  dipole  orientation 
and  altitude  given  by  7  =  17.84°,  4>  =  -58.69°  and 
z  =  12,945  ft. 
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Figure  25.  EY  field  strength  comparisons  between 

'TWIRE"  ( - )  and  point  dipole  (•  •  ■  •)  results 

for  an  azimuth  of  270°.  Point  dipole  orientation 
and  altitude  given  by  'y  =  70.33°,  ^  =  7.23"  and 
z  =  17,250  ft. 
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Figure  26.  EZ  field  strengths  calculated  using 

"TWIRE"  for  azimuths  of  90°  ( - )  and 

270°  ( - ). 


Figure  28.  Rotational  dependence  of  EZ  and  EY  field  strengths  calculated 
using  ’TWIRE". 
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Figure  29.  EZ  field  strength  comparisons  between 

'TWIRE"  ( - )  and  point  dipole  (----)  results 

for  a  laterally  inhomogeneous  guide. 
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Figure  30.  EY  field  strength  comparisons  between 

'TWIRE”  ( - )  and  point  dipole  ( - )  results 

for  a  laterally  inhomogeneous  guide. 
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APPENDIX  A-PROGRAM  LISTING  OF  “TWIRE” 


A-1 


c  twire;  a  program  assembled  from  fas tmc( Ferguson  and  Snyder-NOSC  ID 

c  400, Nov.  1980), from  a  towline  configuration  program (Huang -NADC  AM 
c  6849, June  1969)  and  from  a  radiation  resistance  program(Pappert- 
c  NOSC  TR  1112, June  1986) 
c 

parameter  (maxmds=30 ,nstore=(maxmds*(maxmds+3) )*2+4) 
parameter  (maxsgs=50) 
c 


implicit  real*8  (a-h,o-z) 

common/mc  inpt/xtra(3 ,maxmds) , temp(maxmds) , freq , omega, waveno , mik, 

$  alpha , ak , akl3 , kal3 , ka23 , ngsq , nthsq , sumO , aconst , econst 

common/mc  stor/a(maxmds ,maxmds) , fofr (maxmds) , s(maxmds) , tp(maxmds) , 
$  xval,sigma,epsr ,nrmode,nrslab 

common/mc  plot/dmin, dmax, dtic ,ramin, ramax , ratic , 

$  ampmin, ampmax,amptic ,phsmin,phsmax,phstic , 

$  plotid,bumpid, iopt, frq, power , rcomp , 

$  incl , theta, dst , 

$  talt , al , aalt , ralt, ident ,nrpts , sizex, sizey , cdate(3) , 

$  c time (2) ,nrcrvs ,nrplts ,nrcurv,naplot ,npplot ,npdif f ,nprint 

common  sumout(501) ,xy(501) ,amp(501) ,phs(501) , 

$  hgt ( 3 , maxmds ) , hgr ( 3 , maxmds ) 

dimension  store(nstore) ,dist(20) , idummy(8) 
complex  *16  solsav 

common/com01/antlen(maxsgs) , antrad(maxsgs) , antalt(maxsgs) , 

$  antang(maxsgs) ,npts 

common/com08/solsav(maxmds ,maxsgs) 
equivalence  (store, a) 
integer  rcomp, totape,ctime,cdate 
character*  4  label, ident 
character*40  fname , prof il , blank 
character*80  bcd,plotid,bumpid 

real*4  xy , amp ,phs ,dmin, dmax, dtic , ramin, ramax, ratic , dist , dst , 

$  ampmin , ampmax , amptic , phsmin , phsmax , phs tic , 

$  sizex , sizey , frq , al , talt , aalt , ralt , incl , theta , 

$  tlong, tlat , rbear , power 

real*8  kal3 ,ka23 , nthsq 
complex*  8  sumout 

complex*16  tp , a, s , fofr ,xtra , temp , hgt , hgr , mi , 

$  mik, ngsq , tmpl , tmp2 , tmp3 , tmp4, ratio(4) ,xtra0 (maxmds) 


c 

c  rcomp 
c  power 
c  talt 
c  al 

c  incl, theta 
c  topht 
c  naplot 

c  ampmin, ampmax, amp inc 
c  npplot 

c  phsmin, phsmax, phs inc 
c  sizex, sizey 
c  nrcurv 
c  nprint 
c  -0 

c  -1 


received  component,  =1  is  z,  =2  is  y 
radiated  power  in  kw 
platform  altitude  in  km 

antenna  length  as  percentage  of  wavelength 
orientation  of  transmitter  antenna 
height  of  the  bottom  of  the  ionosphere  in  km 
amplitude  plot  flag 

amplitude  plot  range  and  tic  interval  in  db 
phase  plot  flag 

phase  plot  range  and  tic  interval  in  degrees 

axis  lengths  in  inches 

number  of  curves  per  graph,  max  of  4 

printout  flag 

none 

amp  and  phs  only 


A-2 


-2 

-3 

-4 


rfacmset  cards  and  mode  constants 

conversion  coefficients  and  relative  excitation 

slab  integrals 


iopt-1 


fields  vs  distance 


ralt 

npdiff 

nrptsl 

drain ,  draax ,  dtic 
totape 
tlong, tlat 
rbear 


receiver  altitude  in  km 

phase  difference  plot  flag 

number  of  points  to  output,  max  of  501 

distance  range  and  tic  interval  in  mm 

output  flag,  —1  writes  to  logical  unit  2 

transmitter  location  in  degrees  west, north 

path  bearing  in  degrees  east  of  north 


note:  tlong, tlat, rbear  are  for  totape=l  output  only 


iopt=2 


fields  vs  altitude 


ramin, ramax , ratic  receiver  altitude  range  in  km 

dist,nrd  receiver  distances  in  mm,  number  of  distances 


namelist/datum/iopt ,rcomp , al , incl , theta, talt, ralt , topht , 

$  naplot .npplot , npdiff .nprint ,nrcurv, sizex, sizey, 

$  ampmax , ampmin , amptic , phsmax , phsmin , phs t ic , 

$  dmax , dmin , dtic , tlong , tlat , rbear , power , totape , nrptsl , 

$  ramin , ramax , ratic , dist , nrd , nrpts2 

data  label/' zyx  '/,nf lag/0/, idummy/8*0/, 

$  epslnO/8 . 85434d- 12/ , dtr/ . 01745329252d0/ , 

$  topht/90 . dO/, 

$  tlong , tlat , rbear/0 . , 0 . , 0 . / . totape/0/ , 

$  dist/20*0./.nrd/0/, 

$  nrptsl/501/,nrpts2/51/, 

$  blank/'  '/ 

data  mi/(0 . dO , -1 . dO)/ 

unit  use 

1  primary  input 

2  totape-1  output 

3  temporary  storage  of  me  data 

4  alternate  input  of  data  for  a  single  radial 


bumpid(l :4)— '  ' 

print  2000 
read  2001 , fname 
lunitl-l 

open(unit-l , f ile-fname , status-' old' ) 
open(unit-3 , status-' scratch' , form-' unformatted' ) 


get  date  and  time  for  plot  identification 
call  date(cdate) 
call  time(ctime) 


c  print  1004 , cdate , ctime 

c 

c  get  control  card 

10  continue 

c  print  2002 

read(lunitl , 1002 , end-998)  bed 
print  1003, bed 
lunit4=lunitl 
c 

if(bcd(l:4)  .eq.  'name'  .or.  bcd(l;4)  .eq.  'NAME')  then 
c  print  2003 

c  read  namelist 

read(lunitl , datum) 

print  datum 

idummy(3)=int(al+. 5) 

if (naplot+npplot+npdiff  .eq.  3)  then 

print  *, 'NAPLOT,NPPLOT,NPDIFF  can  not  all  be  1' 
go  to  10 
endif 

ifCphsmin  .eq.  0.  .or.  phsmax  .eq.  0.)  then 
print  *,'PHSMIN  or  PHSMAX  is  zero' 
go  to  10 
end  if 

lf(nrcurv  .gt.  4)  then 
print  *,'NRCURV  .gt.  4' 
go  to  10 
end  if 

if(sizey  .gt.  9.)  then 
print  *,'SIZEY  .gt.  9' 
go  to  10 
end  if 

if(iopt  .eq.  1)  then 

if(nrptsl  .gt.  501)  then 
print  *,'NRPTS1  .gt,  501' 
go  to  10 
end  if 

nrpts-nrptsl 

else 

nrpts=nrpts2 
end  if 
go  to  10 
end  if 
c 

if(bcd(l:4)  .eq.  'bump'  .or.  bcd(l;4)  .eq.  'BUMP')  then 
c  store  bump  id  card 

bumpid=bcd 
go  to  10 
end  if 
c 

if(bcd(l:4)  .eq.  'outp'  .or.  bcd(l:4)  .eq.  'OUTP')  then 
J»7 

do  while  (bcd(j:J)  .eq.  '  ') 
if(j  .eq.  80)  then 

print  *, 'Output  file  name  not  found' 
go  to  10 


end  if 

j-j+1 

end  do 

fname=bcd( j : ) 

open (uni t=2 , f ile=fnanie , status- 'new' , form- 'unformatted' ) 
go  to  10 
end  if 
c 

if(bcd(l;4)  .eq.  'inpu'  .or.  bcd(l:4)  .eq.  'INPU')  then 
j-7 

do  while  (bcd(i:i)  .eq.  '  ') 
if(j  .eq.  80)  then 

print  *, 'Input  file  name  not  found' 
go  to  10 
end  if 

j-j+1 

end  do 

fname=bcd( j : ) 

open(unit— 4 , f ile— fname , status—' old' ) 
lunit4=4 

read(lunit4 , 1002)  bed 
print  1003, bed 
go  to  19 
end  if 
e 

if(bed(l;4)  .eq.  'elos'  .or.  bed(l:4)  .eq.  'CLOS')  then 
if(bod(7:7)  .eq.  'o'  .or.  bed(7;7)  .eq.  '0')  then 
olose(unit-2) 
else 

elose (unit-4) 
lunit4=lunitl 
end  if 
go  to  10 
end  if 
e 

19  if(bcd(l;4)  .eq.  'sw  '  .or.  bod(l:4)  .eq.  'SW  ')  then 
e  get  xmtr  and  profile  parameters 

read(bod, 1018)  tlong,tlat,rbear,profil 
go  to  20 
end  if 
e 

if(bed(l;4)  .eq.  'ipsq'  .or.  bed(l:4)  .eq.  'IPSQ')  then 
e  get  xmtr  and  profile  parameters 

read (bed, 1018)  tlong, tlat , rbear ,prof il 
go  to  20 
end  if 
c 

lf(bcd(l:4)  .eq.  'data'  .or.  bcd(l:4)  .eq.  'DATA')  then 
prof il-blank 
rbear-720 . 
go  to  20 
end  if 
c 

if(bcd(l:4)  .eq.  'star'  .or.  bcd(l:4)  .eq.  'STAR')  go  to  30 
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if(bcd(l:4)  .eq.  'quit'  .or.  bcd(l:4)  .eq.  'QUIT')  go  to  998 
c 

print  *, 'Error  in  control  card' 
go  to  10 
c 

c  read  propagation  path  data 

20  read(lunit4 , 1002)  plotid 

if(rbear  .ne.  720.)  then 

do  while  (plotid(j : j+12)  .ne.  '  ') 

if(j  .eq.  67)  then 

print  *,'Plot  identification  is  too  long  to  append  bearing' 
go  to  22 
end  if 

j-j+1 

end  do 

write(plotid(j : j+12) , 1021)  rbear 
end  if 

22  print  1003, plotid 
nrmode-maxmds 
nrslab=0 
rho=-l. 

th==0 . 

nthsq=l . +alpha*topht 
rewind  3 

23  read(lunit4 , 1002)  bed 

read(bcd , 1020)  rr , f f , aa , cc , bb , ss , ee 
if(rr  .ne.  40.  .and.  ss  .eq.  0.)  go  to  23 
xval-rr*1000 . 

if(nrslab  .gt.  0)  write(3)  store 
c 

if(rr  .It.  40.)  then 
nrslab=nrslab+l 
nm-O 

bb-bb*10000 . 

if(bcd(71;71)  .eq.  't'  .or.  bcd(71:71)  .eq.  'T') 

$  read(bcd, 1028)  th 

if(nprint  .ge.  1)  print  1022 ,nrslab , rr , f f , aa , cc ,bb , ss , ee , th 
if(nprint  .ge.  4)  print  1024 
if(rr  .eq.  0.)  then 

c  if(nprint  . le .  1)  print  1022 .nrslab , rr , ff , aa , cc ,bb , ss , ee , th 

frq=ff 
freq=f f 

const-682 . 2408*dsqrt(freq) 
omega-6 . 283185308d3*freq 
waveno=20 . 958445d-3*freq 
mik-dcmplx(0 . dO , -waveno) 
aconst--8 . 686d3*waveno 
econs t=20 . *dlogl0 (35 .*waveno) 
ak=a 1 pha/waveno 
akl3-=dexp  (dlog(ak)/3 .  dO) 
kal3-l.d0/akl3 
ka23=kal3**2 
end  if 

if(rho  .gt.  rr)  then 
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print  *,'Rhos  are  out  of  order' 
go  to  10 
end  if 
rho-rr 
sigraa“ss 
epsr-ee 

ngsq“dcmplx(epsr , -sigma/(omega*epslnO) ) 
if(th  .gt.  0.)  nthsq-1 .+alpha*th 
24  read(lunit4, 1023)  indxl , trl , til , itrml , tmpl , tmp2 

c 

if (indxl  .gt.  0)  then 

read(lunit4 , 1023)  indx2 , tr2 , ti2 , itrm2 , tmp3 , tmp4 
if(nin  .eq.  maxmds)  go  to  24 
if (zabs(tmpl)  .eq.  0.)  go  to  24 

if(trl  .ne.  tr2  .or.  til  .ne.  ti2  .or.  itrml  .ne.  itrm2)  then 
print  *, 'Modes  are  out  of  order' 
go  to  10 
end  if 

if(itrml  .eq.  0)  then 

print  *, 'Mode  conversion  flag  is  missing' 
go  to  10 
end  if 
nm— nm+1 

if(nprint  .ge.  4)  then 

print  1025 ,nm, indxl , trl , til , itrml , tmpl , tmp2 , 

$  indx2 , tr2 , ti2 , itrm2 , tmp3 , tmp4 

end  if 

tp(nm)-dcmplx(trl ,til) 
s(nm)“zsln(tp(nm)*dtr) 
rat io(2*indxl-l) -tmpl 
ratio(2*indxl  )-tmp2 
ratio (2*indx2-l)-tmp3 
ratio(2*indx2  )-tmp4 
temp (nm) -ratio (4) 

c  get  vertical  excitation  factor  for  printout 

xtra0(nra)-( -2 . 124292961d0,0 . d0)*ratio(l)*s (nm)**2 
c  get  ey/hy 

if (itrml  .eq.  1)  then 

fofr (nra)-ratio(3)/ratio(l) 
else 

fofr (nm)-ratio(2)/(ratio(3)*ratio(4) ) 
end  if 

if(nrslab  .eq.  1)  then 

c  get  hy  excitation  factor  at  the  transmitter 

xtra(l ,nm)-- ratio (l)*s(nra) 
xtra(2,nra)=  ratio(3)*ratio(4) 
xtra(3,nm)=  ratio(l) 
end  if 
go  to  24 
else 

if(nprint  .gt.  1)  print  1027, nm 
nrmode-nra 
end  if 
c 

if(nprint  .gt.  1)  then 
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c 


26 

c 

c 


c 

c 


c 

c 


30 


c 


c 


31 

c 

c 


c 

c 

c 

c 


print  mode  constants  table 
print  1031 
do  26  m-l.nrmode 
sr=dble(s(m) ) 
si-dimag(s (m) ) 
atten=aconst*si 
voverc=l . dO/sr 
wr=dble (xtraO (m) ) 
wi=diraag(xtraO (m) ) 

wml=econs t+10 . d0*dlogl0(dmaxl (wr**2+wi**2 , 1 . d-20) ) 
wal=datan2 (wi , wr) -1 . 5707963267949d0 
print  1032 ,m, atten, voverc , wml , wal 
continue 
end  if 

get  conversion  coefficients 
call  mcstep 
go  to  23 

else 

if(nrslab  . le .  1)  then 

print  Insufficient  number  of  slabs, 

$  'horizontally  homogeneous  only' 

end  if 
end  if 

calculate  mode  sum 

ident”label(rcomp:rcomp)//'  ' 

sumO=const*sqrt (power) 

call  tacamo 

call  coord 

call  rpower 

if(iopt  .eq.  1)  then 

deltad=(dmax-dmin)/float(nrpts-l) 

dst=dmin 

do  31  i-l,nrpts 

xy(i)~dst 

call  me fids ( i , dst , sumout(i) , amp( i) ,phs ( i) ) 
ds  t=ds  t+de 1 tad 

if(totape  .gt.  0) 
save  mode  sum 

$  write(2)  sumout , frq , tlong, tlat, rbear , power , rcomp , incl , theta , 
$  talt,ralt,dmin,dmax, idummy.profil 

mode  sum  output 
call  mcplts 

else 

deltar<=(ramax-ramin)/float(nrpts-l) 
do  43  i=l,nrd 
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41 


42 

43 


49 


998 


1000 

1002 

1003 

cl004 

1018 

1020 

1021 

1022 


dst-dist( i) 
ralt-0 . 

call  mcflds (i , dst ,  suniout(l)  , amp(l)  ,phs(l)  ) 
phsl=phs (1) 
cycle=0 . 

ralt=ramin 
do  42  j“l,nrpts 

call  htgain(2 , freq , sigma, epsr , alpha, nrmode , tp , dble(ralt) ,hgr) 
tnipl-(0.d0,0.d0) 
do  41  m=l, nrmode 
do  41  jj=l,npts 

tmpl=tmpl+solsav(m, j j )*hgr(rcomp ,m) 

tar=tmpl 

tai=tmpl*mi 

phs2-=datan2  (tai  ,  tar)*57 . 2957795d0 
if (dabs(phsl-phs2)  .ge.  180. dO)  then 
if(phsl  . le .  phs2)  then 
cycle-cycle -  360 . dO 
else 

cycle-cycle+360 . dO 
end  if 
end  if 
phsl=phs2 
xy(j )=ralt 

amp( j )-10 .*dlogl0( tar**2+tai**2) 

phs ( j )=phs2+cycle 

sumout(j  )*=tmpl 

ralt=ralt+deltar 

call  mcplts 

ralt-0. 

call  htgain(l , freq , sigma, epsr , alpha, nrmode , tp , dble(ralt) ,hgr) 
end  if 

if (nprint.gt.O)  print  1000 
if(lunit4  .eq.  4)  then 

read(lunit4 , 1002 , end-49)  bed 
print  1003, bed 
go  to  19 
close(unit-4) 
end  if 
go  to  10 

print  1998 , nrcrvs , nrplts 
stop 

format ( ' 1 ' ) 
format(a80) 
format(lx,a80) 

formate 'Additional  plot  identification:  ’ , 3a4 , lx , 2a4) 
format (9x , f 7 . 0 , 2f6 . 0 , 6x, a40) 

format(lx,f7.0,3(2x,f8.0),2(2x,el0.0) ,2(2x,e5.0)) 
formate ' .  rbear-' , f5 . 1) 

formate'  Slab  ',i2,'  R',f7.3,'  F',f8.4,'  A',f8.3,'  C',f8.3,'  M' , 
1  f6.3,'  S',lpel0.3,'  E',0pf5.1,'  T',f5.1) 
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1023 

1024 

1025 

1026 

1027 

1028 

1031 

1032 
1998 
2000 
2001 
2002 
2003 


$ 

$ 


$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 

$ 


format (il,2f 9. 0,il,4el5.0) 

format(/llx,  '  M  ID  THETA'") 

format(llx,i2,3x,il,0p2fl0.5,i2,2(lx,lp2el6.8)/ 

16x,  il,0p2fl0.5,i2,2(lx,lp2el6.8)) 

format (16x,  il ,0p2fl0. 5 , 12 , 2(lx, Ip2el6 .8)/ 

16x,  il,0p2fl0.5,i2,2(lx,lp2el6.8)) 

format( ' , 80x , '  Modes ’,13) 
format (7 lx, f5 . 0) 

formate'  MODE  ATTEN  V/C  WAIT"S  EXC') 

format (14 , f8 . 3 , f9 . 5 , 2 (f 10 . 3 , f7 . 3) ) 

formate '  End  of  job.  ',12,'  curves  plotted  in  ',12,'  graphs') 
formate'  Enter  .com  file  name:  ') 
format (a40) 

formate'  Control  card:  [name  input  output  start  quit  close]  '$) 
formate'  Enter  DATUM:  iopt[l]  rcomp[l]  power [1]'/ 

'  incl[0]  theta[0]  talt[0]  ralt[0] '/ 

'  dmin[0]  dmax[10]  dtic[l]'/ 

'  tlong[0]  tlat[0]  rbear[0]  totape[0]'/ 

'  ampmin[-50]  ampmax[70]  amptic[10]'/ 

’  phsmin[-360]  phsmax[360]  phstlc[90]'/ 

'  si2ex[5]  sizey[6]'/ 

'  nrd[0]  dist[20*0]'/ 

'  ramin[0]  ramax[50]  ratic[5]'/ 

'  nrpts[501]  nrptsl[501]  nrpts2[51]'/ 

'  nrcurv[l]  naplot[l]  npplot[0]  nprint[0]') 


end 
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</></>■</> -C/>  </></></></> -Crt  </>  </>  </></></>  </>-(/></></>-</> 


block  data 


parameter  (maxmds-30 ,nstore-(maxmds*(maxmds+3) )*2+4) 
implicit  real*8  (a-h,o-z) 

common/mc  inpt/xtra(3 .maxmds) , temp(maxmds) ,freq, omega, waveno.mik, 
$  alpha, ak, akl3 ,kal3 ,ka23 ,ngsq,iithsq , sumO , aeons t , econst 

common/mc  plo t/dmin , dmax , dtic , ramin , ramax , ratic , 
ampmin , ampmax , amp tic , phsmin , phsmax , phs tic , 
plotid,bumpid, iopt, frq, power , rcomp , 
incl , theta, dst , 

talt , al , aalt ,ralt , ident .nrpts , sizex, sizey , cdate(3) , 
ct ime (2 ) , nrcrvs , nrplts , nrcurv, naplot , npplot , npdif f , nprint 
integer  rcomp, ctime,cdate 
character*  4  ident 
character*80  plotid,bumpid 
real*8  kal3 ,ka23 ,nthsq 
complex*16  xtra, temp,mik,ngsq 
real*4  dmin, dmax, dtic , ramin, ramax, ratic , dst , 

ampmin , ampmax , amptic , phsmin , phsmax , phstic , 
sizex, sizey , frq , al , talt, aalt , ralt , incl , theta, 
power 

data  iopt/1/, 

alpha/3 . 14d-04/ , 

rcomp/1/, incl, theta/0. ,0./. talt,ralt/0. ,0./. 

power/1 . / , al/0 . / , 

dmin, dmax, dtic/0. ,10. ,1./, 

ramin , ramax , ratic/0 . , 50 . , 5 . / , 

ampmin, ampmax, amp tic/ -50. ,70. ,10./, 

phsmin, phsmax, phstic/- 360. ,360. ,90./, 

sizex, sizey/5 . ,6./, 

nrpts/501/, 

nrcurv/1/ , naplo t/1/ , npplot/0/ , npdif f/0/ , nprint/0/ , 
nrcrvs ,nrplts/0 ,0/ 

end 
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subroutine  coord 


implicit  real*8  (a-h,o-z) 
parameter  (maxsgs-50) 
real*4  gamma, phi 

common/com01/antlen(maxsgs) , antrad(maxsgs) , antalt(maxsgs) , 

$  antang(maxsgs) ,npts 

common/com02/cpl ,psi , iclock, chicp 
common/com03/x(maxsgs) ,y(maxsgs) ,z(maxsgs) 
common/com04/xmid(maxsgs) ,)nnid(maxsgs) ,zmid(maxsgs) 
common/com05/gamma (maxsgs ) , phi (maxsgs ) 

nsegs=npts-l 

calculate  coordinates  of  ends  of  segments  of  the  antenna 
format(2x, '  segment' ,9x, 'x' ,llx, 'y' ,llx, 'z' ) 
do  i=l,nsegs+l 

X ( i ) =antr ad ( i ) *co  s ( antang ( i ) ) 
y(i)=antrad(i)*sin(antang(i) ) 
z(i)=antalt(i) 
enddo 

Transformation 

capphi=datan2 (y(npts) ,x(npts)) 
do  i<=l,nsegs+l 

xprime-  x(i)*cos(capphi)+y(i)*sin(capphi) 
yprime--x(i)*sin(capphi)+y(i)*cos (capphi) 
x( i)-xprime 
y(i)-yprime 
enddo 

do  i-l,nsegs+l 

x(i)«-x(i)-x(npts) 
y(i)=y(i) -y(npts) 
enddo 

if  airplane  is  moving  in  a  clockwise  orbit  (iclock=l)  then  change  y 
if (iclock  .eq.  1)  then 
do  i-l ,ns '"s+l 

y(i)-=-y(i) 

enddo 

endif 

cpsi=cos(psi*0 .0174533) 
spsi=sin(psi*0 . 0174533) 
do  i=-l,nsegs+l 

xpp”  x(i)*cpsi+y(i)*spsi 
ypp=-x(i)*spsi+y(i)*cpsi 
x(i)”xpp 
y(i)-ypp 
enddo 

z(nsegs+2)-z(nsegs+l) -cpl*sin( chicp* .0174533) 


x(nsegs+2)-  -cpl*spsi*cos(chicp* . 0174533) 
y(nsegs+2)--cpl*cpsi*cos (chicp* . 0174533) 
if(iclock  .eq.  1)  then 
x(nsegs+2)— -x(nsegs+2) 
y(nsegs+2)--y(nsegs+2) 
endif 

do  i“l,nsegs+l 

if(antlen(i)  .le,  (antlen(nsegs+l)+cpl)/2 .  .and. 

$  antlen(i+l)  .gt.  (antlen(nsegs+l)+cpl)/2 . ) then 

iniid=i 
go  to  915 
end  if 
end  do 

915  print  *,  ' iraid=' , imid 

c  calculate  midpoints  of  segments  of  the  antenna 

do  i-l,nsegs+l 

xmid(i)-(x(i+l)+x(i) )/2 .0 
ymid(i)-(y(i+l)+y(i))/2.0 
zmid(i)“(z(i+l)+z(i))/2.0 
end  do 

c  calculate  gamma  and  phi 

c  gamma=inclination  angle  of  antenna  from  the  vertical 

c  phi=azimuth 

print  AFTER  TRANSFORMATION' 
print  910 

910  format(2x, '  segment' ,6x, 'xmid' ,8x, 'ymid' ,8x, 'zmid' , lOx, 'gamma' , 
$  7x,'phi') 

do  i-l,nsegs+l 
delx-x(i+l) -x(i) 
dely-y(i+l)-y(i) 
delz-z(i+l)-z(i) 

geimma  (  i )  =atan2  ( s  qr  t  (  de  lx**2+de  ly**2  )  ,  de  1  z ) 

phi(i)“atan2(dely ,delx) 

print  905 , i ,xmid(i) ,ymid(i) ,zmid(i) , 

$  gamma(i)*57. 29578, phi(i)*57. 29578 

905  format('  ' ,4x, i3 , 6x, 5(fl0. 3 , 2x) ) 
enddo 
return 
end 
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subroutine  mcf Ids ( index , ds t , sumout , amp , phs ) 
c 

c  calculate  mode  sum 
c 

parameter  (maxmds-30 ,nstore-(maxmds*(maxmds+3) )*2+4) 
parameter  (maxsgs=50) 
c 

implicit  real*8  (a-h,o-z) 
real*8  idl 
real*4  gamma, phi 
c 

common/com01/antlen(maxsgs) , antrad(maxsgs) , antalt(maxsgs) , 

$  antang(maxsgs) ,npts 

coramon/com04/xmid(maxsgs) ,ymid(maxsgs) , zmid(maxsgs) 
common/com05/gamma(maxsgs) ,phi(maxsgs) 
common/com06/idl (maxsgs) 

common/mc  inpt/xtra(3 .maxmds) , temp(maxmds) , freq , omega .waveno ,mik, 

$  alpha, skip (7) ,constO , aconst , econst 

common/mc  stor/a(maxmds , maxmds) , fofr (maxmds) ,s(maxmds) ,tp(maxmds) , 
$  xtwo , sigma, epsr ,nrmode ,nrslab 

common/mc  plot/iskipl (12) ,plotid(20) , iskip2 (23) , rcomp , 

$  incl , theta , iskip3 , talt , al , aalt , ralt , 

$  ident , iskip4(14) ,nprint 

common  iskip5(2505) ,hgt(3 ,maxmds) , hgr ( 3 , maxmds ) , 

$  soln  a(maxmds , maxsgs) , soln  b (maxmds , maxsgs) 

dimension  store(nstore) , t(3) 
complex*16  ssav 
dimension  ssav(maxmds) 
equivalence  (store, a) 

complex*16  a , s , tp , fofr , xtra , temp , solna , solnb ,hgt ,hgr , 

$  ta, tb ,mi ,one ,mik,mikx 

complex  *16  solsav 
common/com08/sol3av(maxmds , maxsgs ) 
real*4  amp , phs , incl , theta , talt , al , aalt , ralt , dst 
integer  rcomp , tcomp .pnrnode 
complex*  8  sumout 
complex*16  sum 
character*4  plotid, ident 
data  mi/(0.d0, -I.d0)/.one/(l.d0,0.d0)/ 
c 

constl=2 . 853d-3*freq**l . 5 
sum— 0 . 

rho=dst*1000 . 
numseg-npts - 1 

900  format('  i, dst, ant  len, rad, alt , ang=' , i3 , 5el2 . 5) 

905  format('  i , dst ,x,y , zmid, gamma.phi , idl=' , i3 , 7el2 . 5) 
c 

if (index  .eq.  1)  then 
rewind  3 
read(3)  store 
do  m-1 , nrmode 
ssav(m)-s (m) 
enddo 
xone-0 . dO 
phs 1-0 . 
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cycle-0 . 
mikx-mik*xtwo 

call  htgain(l , freq , sigma, epsr, alpha, nrmode , tp , dble(ralt) ,hgr) 
end  if 

do  i-1 , numseg+1 

if (index  .eq.  1)  then 
c 

c  antenna  orientation  factors 

t(l)-cos(gamma(i) ) 
t(2)-sin(gamma(i) )*sin(phi(i) ) 
t(3)-sin(ganima(i)  )*cos(phi(i) ) 
c 

c  xone  is  distance  of  previous  slab  from  the  transmitter 

c  xtwo  is  distance  of  next  slab 

c 

c  get  data  for  first  slab 

call  htgain(l, freq, sigma, epsr, alpha, nrmode, tp,antalt(i) ,hgt) 
do  689  m-1, nrmode 
c  hy  excitation  factor 

ta-(0.d0,0.d0) 
do  686  tcomp-1,3 

686  ta-ta+xtra(tcomp ,m)*hgt(tcomp ,m)*t(tcomp) 

soln  a(m, i)-ta 
if(rcomp  .eq.  1)  then 
c  ez  excitation  factor 

ta“-s(m)*ta 
else 

c  ey  excitation  factor 

ta-fofr{m)*ta 
end  if 

soln  b(m, i)-ta*hgr (rcomp ,m) 

689  continue 

c 

end  if 
enddo 
c 

if(dst  .eq.  0.)  then 
amp-100 . dO 
phs-0 . 
sumout-amp 
return 
end  if 
c 

720  if(rh_  "t.  y:*'wo)  then 
c 

c  end  of  current  slab 

mikx-mik* (xtwo - xone ) 
pnmode— nrmode 
do  i-1, numseg+1 
do  735  m=l, pnmode 

c  get  excitation  factors  at  end  of  slab 

soln  a(m,i)-soln  a(ra, i)*zexp(mikx*(s (m) -one) ) 

735  solsav(m, i)-soln  a(ra,i) 

enddo 
xone-xtwo 
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c 

c  get  data  for  next  slab 

read(3)  store 
mikx— mik*xtwo 

call  htgain(l , freq , sigma ,epsr, alpha, nrmode , tp , dble(ralt) ,hgr) 
do  i— l,numseg+l 

do  740  m2-l, nrmode 
c  hy  excitation  factor 

ta=(0 . dO , 0 . dO) 
do  737  ml=l,pnmode 
737  ta=ta+solsav(ml , i)*a(ml ,m2) 

soln  a(m2,i)=ta 
if(rcomp  .eq.  1)  then 
c  ez  excitation  factor 

ta=-s(m2)*ta 
else 

c  ey  excitation  factor 

ta=fofr(m2)*ta 
end  if 

soln  b(m2 , l)-ta*hgr(rcomp,m2) 

740  continue 

c 

enddo 
go  to  720 
c 

else 

c 

c  get  sum  of  modes 

mikx"=mik*  ( rho  -  xone  ) 
do  i"=l  ,numseg+l 
ta=0 . dO 

factor=*idl(i)*constl/dsqrt(dabs(dsin(  (rho-xniid(i)  )/6 . 366d3) ) ) 
do  730  m-1, nrmode 

tb=solnb(m, i)*zexp(mikx*(s(m) -one) )*zexp( -mik*xmid(i)*ssav(m) ) 
$  *factor 

solsav(m, i)=tb 
730  ta-ta+tb 

sum-sum+ta 
enddo 
end  if 
sumr=sura 
sumi=mi*sum 

phs2-datan2 (sumi , sumr)*57 . 2957795d0 
if (dabs(phsl-phs2)  .ge.  180. dO)  then 
if(phsl  .le.  phs2)  then 
cycle-cycle -  360 . dO 
else 

cycle=cycle+360 . dO 
end  if 
end  if 
phsl-phs2 
sumout-sum 

amp-10 .*dlogl0(dmaxl(l . d-30 , sumr**2+sumi**2) ) 

phs-phs2+cycle 

return 
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c 

c 

1030  format( ' 1' , 20a4/lx,al , '  component  INCL  -',f4.0, 

$  '  THETA  ,f 5.0, '  PALT  '  AALT-',f5.1. 

$  '  RALT  .f5.1) 

1031  formate 'OSLAB  RHO  MODE  ATTEN  V/C  REL  EXC  1'/ 

$  Ix,i3,f9.3) 

1032  format ( 13x , 13 , f 8 . 3 . f 9 . 5 . 2 (f 10 . 3 , f 7 . 3 ) ) 
end 
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</>  </>  </>  </>  </>  ^  ^  ^ 


c 

c 

c 


c 

c 


c 


19 

c 

c 


subroutine  mcplts 

plot  mode  sum  amplitude  and  relative  phase 


real*8  cpl ,psi , v, rpl , zpl , zplane , radius , chicp 
common/com02/cpl , psi , iclock, chicp 
common/com09/v, rpl , zpl 

common/mc  plo t/dmin , dmax , dtic , ramin , ramax , ratio , 
ampmin , ampmax , amp  tic, phsmin , phsmax , phs  t i c , 
id(20) ,bumpid(20) , iopt.freq, power, rcomp, 
angl , ang2 , dst , talt , 

al , aalt , ralt , ident ,nrpts , sizex, sizey , cdate(3) , c time (2) , 
nrcrvs ,nrplts ,nrcurv,naplot ,npplot ,npdiff ,nprint 
common  skipl(1002) ,xy(501) ,outl(501) ,out2(501) ,phil(501) ,up(501) 
dimension  xl(2) ,yl(2) ,ul(2) . label(7 , 5) , iplbll(13) , iplbl2(13) 
integer  bumpid, cdate , ctime , rcomp ,headng( 13) 
logical  up ,ul ,nframe ,pflag 
character*  4  ident 
character*28  hlabel(5) 
character*52  hhead, blank, p Ibl 1 ,plb 12 
character*80  plotid, saveid 
character*80  bumpch 

equivalence  ( label ,hlabel) , (headng, hhead) , (id, plotid) , 
(iplbll,plbll) ,(iplbl2,plbl2) 
equivalence  (bumpch, bumpid) 
data  nframe/. true ./ 
data  jopt/0/ 

data  xl/'l.l, -0.1/, ul/2*. false./ 

data  hhead/'  Freq  Pwr  Aalt  Ralt  Rpl  V  Psi  Iclock  '/ 
data  blank/'  '/ 

data  hlabel/ 

'  Amplitude  (dB  above  luv/m)  ' , '  Relative  phase  (Degrees)  ' , 

'  Phase-phasel  (Degrees)  ','  Distance  (Megameters)  ', 

'  Altitude  (kilometers)  '/ 


convert  plane  altitude  and  plane  radius  from  ft  to  km 
zplane=zpl*l . Od-3 
radius-rpl*! . Od- 3 

convert  receiver  altitude  from  km  to  kft 
ral-ralt/0.3048 
if(nprint  .gt.  0)  then 
if(iopt  .eq.  1)  then 
print  1001 , ident , ralt 
else 

print  1003 , ident , dst 
end  if 

if(outl(l)  .ge.  100.)  outl(l)-99.99 

nl-nrpts/4+1 

do  19  il=l,nl 

print  1011 , (xy( i) , outl (i) ,out2(i) , i=il ,nrpts ,nl) 
end  if 


if (naplot+npplot+npdiff  .eq.  0)  return 


if (j opt  .ne.  iopt)  then 
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jopt-iopt 
ndata-0 
line-0 
end  if 

ndata-ndata+1 
line-line+1 
if(ndata  .eq.  1)  then 
if (.not.  nframe)  then 
call  pltend 
line— 0 
end  if 

nframe- . true . 
do  1  i-l,nrpts 
1  phil(i)-out2(i) 

end  if 

if(iopt  .eq.  1)  then 
xmin-dmin 
xmax—draax 
xtic-dtic 
else 

ymin-ramin 
ymax-ramax 
ytic-ratic 
end  if 
c 

if(naplot  .eq.  1)  then 
nplot-1 
amax-ampmax 
arain-ampmin 
atic-amptlc 
do  29  i-l,nrpts 
if(outl(i)  .le.  amax)  then 
ff(outl(i)  .ge.  amin)  then 
up(i)-. false, 
else 

up  (i)-. true . 
outl(i)=amin 
end  if 
else 

up  (i)-. true . 
outl(i)-amax 
end  if 

29  continue 

if(iopt  .eq.  1)  then 
ymax—araax 
ymin-amin 
ytic-atic 
else 

xmax— amax 
xmin— amin 
xtic-atic 
end  if 
go  to  100 
end  if 
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30 


if(npplot  .eq.  1)  then 
np lot-2 
fk-0. 

do  39  i-l,nrpts 
outl(i)-out2(i)+fk 
up(i)-. false . 

35  if(outl(i)  .ge.  phsmin)  then 

if(outl(i)  .le.  phsmax)  go  to  39 
fk=fk-phsmax 
outl(i)— outl(i) -phsmax 
else 

fk-fk-phsmin 
outl(i)-outl(i) -phsmin 
end  if 
up ( i ) = . true , 
go  to  35 

39  continue 

if(iopt  .eq.  1)  then 
ymax=phsmax 
ymin-phsmin 
ytic-phstic 
else 

xmax-phsmax 
xmin— phsmin 
xtic-phstic 
end  if 
go  to  100 
end  if 

c 

40  if(npdiff  .eq.  1  .and.  iopt  .eq.  1  .and.  ndata  .gt.  1)  then 

nplot-3 

fk-0. 

do  49  i-l,nrpts 
outl(i)-out2 (i) -phil(i)+fk 
up( i)-. false . 

45  if(outl(i)  .ge.  phsmin)  then 

if(outl(i)  .le.  phsmax)  go  to  49 
fk-fk- phsmax 
outl(i)=outl(i) -phsmax 
else 

fk-fk -phsmin 
outl(i)-outl(i) -phsmin 
end  if 
up(i)-. true . 
go  to  45 

49  continue 

if(iopt  .eq.  1)  then 
ymax-phsmax 
ymin-phsmin 
ytic-phstic 
else 

xmax-phsmax 
xmin-phsmin 
xtic-phstic 
end  if 
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go  to  100 
end  if 
c 

50  if (line  .eq.  nrcurv)  then 
call  pltend 
nframe-. true . 
line-0 
else 

call  pltoff 
saveid-plotid 
yp0=yp 
end  if 
return 
c 

100  ipwr=power 
iangl-angl 
iang2“ang2 

if(line  .eq.  1)  then 
if (nframe)  then 
call  pltbgn 
nrplts-nrplts+l 
xp-0 . 
yp-0. 

c  call  symbol(xp,yp, .1. 'FASTMC' ,0. ,6) 

xp-xp+0 . 7 

call  number (xp ,yp , . 1 , float (nrplts) , 0 . , -1) 
xp=xp+l . 0 

call  symbol (xp, yp, .l,cdate,0. ,9) 
xp-xp+1 . 0 

call  symbol  (xp,  jrp,  .l,ctime,0.  ,8) 
call  plot(l .  1 ,0. 8+0.45*nrcuirv, -3) 
nframe-. false . 
end  if 

call  bordr2 (sizex,xmin,xmax,xtic , 2 .*xtic , -1 , 

$  sizey ,ymin,ymax,ytic , 2 .*ytic , -1) 

if(iopt  .eq.  1)  then 
xp-- . 7 

yp-  . 5*(sizey-2 . 8) 

call  symbol (xp ,yp , . 1 , label(l ,nplot) ,90. ,28) 
xp-  . 5*(sizex-2 . 8) 
yp—  .4 

call  symbol (xp, yp, .1, label (  1,4  ),  0.,28) 

else 
xp=- . 7 

yp-  . 5*(sizey-2 . 8) 

call  symbol (xp ,yp ,. 1 , label (  1,5  ),90.,28) 

xp-  . 5*(sizex-2 . 8) 
yp—  .4 

call  symbol (xp ,yp ,. 1 , label(l ,nplot) ,  0.,28) 
end  if 
xp-0. 
yp-- .6 

if (bumpch(l :4)  .ne.  '  ')  then 

call  symbol (xp ,yp , .l,bumpid(l) ,0. ,4) 
xp-xp+0 . 5 
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call  symbol (xp , yp , . 1 ,bumpid(2) , 0 . ,68) 
xp“0 . 
yp=yp- .15 
end  if 

call  symbol(xp ,yp , . 1 , id,0 . , 68) 

saveid=plotid 

yp=yp- .15 

call  symbol (xp ,yp 1 .headng, 0 ., 52) 
yp=yp- .15 
dummy=9999 . 

write (plb 11 , 2000) ident , f req , ipwr , zplane , ral , radius , v , psi , iclock 
2000  format (al , lx , f 6 . 2 , lx, i4 , 3(lx, f 5 . 2) , lx, f 5 . 1 , lx , f5 . 1 , lx , i4) 
call  symbol (xp ,yp , . 1 , iplbll , 0 . ,52) 
xp=0 . 
yp0=yp 
else 

call  plton 
yp=yp0- .15 

write (plbl2 , 2000) ident , freq, ipwr , zplane , ral , radius , v, psi , iclock 
c 

c  check  for  common  data  in  plot  label 

if(plbl2(  1:1  )  .eq.  plbll(  1:1  )) 

$  plbl2(  1:1  )='  ' 

if(plbl2(  3:8  )  .eq.  plbll(  3:8  )) 

$  plbl2(  3:8  )=' 

if(plbl2(10:13)  .eq.  plbll(10 : 13) ) 

$  plbl2(10:13)=' 

if(plbl2(15:19)  .eq.  plbll(15 ; 19) ) 

$  plbl2(15:19)=' 

if(plbl2(21:25)  .eq.  plbll(21 : 25) ) 

$  plbl2(21;25)=' 

if(plbl2(27:31)  .eq.  plbll(27 : 31) ) 

$  plbl2(27:31)-' 

if(plbl2(33:37)  .eq.  plbll(33 : 37) ) 

$  plbl2(33:37)='  ' 

if(plbl2(39:43)  .eq.  plbll(39:43)) 

$  plbl2(39:43)=’ 

if(plbl2(45:48)  ,eq.  plbll(45 :48) ) 

$  plbl2(45:48)=' 

c 

if(plotid  .eq.  saveid)  then 
pf lag= . false . 
else 

call  symbol(xp ,yp , . 1 , id, 0 . , 68) 
pflag” . true . 
end  if 

if(plbl2  .ne.  plbll  .and. 

$  plbl2  .ne.  blank)  then 

if(pflag)  then 
yp=yp- . 15 

call  symbol (xp ,yp ,. 1 .headng, 0 52) 
yp-yp- .15 
end  if 

call  symbol (xp ,yp , . 1 , iplbl2 , 0 . ,52) 
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end  if 


yl(l)-yp+.05 

yl(2)=yp+.05 

call  curve(xl ,yl ,ul , 2 ,0.,0.,1.,1. .line) 
c 

c  draw  curves 

nrcrvs=nrcrvs+l 
if(iopt  .eq.  1)  then 

call  curve(xy , outl.up .nrpts .xmin.ymin, (xmax-xmin) /sizex 
$  (yniax-yinin)/sizey ,  line) 

else 

call  curve(outl,xy, up, nrpts, xnjin.ymin,  (xniax-xniin)/sizex 
$  (yinax-yTnin)/sizey ,  line) 

end  if 
c 

if(nplot-2)  30,40,50 
c 

1001  formate'  ',al,'  COMP  RALT  •=  ',f7.3/ 

$  4('  DIST  AMPLITUDE  PHASE', 7x)) 

1003  formate'  '.al,'  COMP  DIST-  '.f6.3/ 

$  4e'  RALT  AMPLITUDE  PHASE', 7x)) 

1011  formate4ef7.3,2fl0.4,5x)) 

end 
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-C/>  </></></>-(/>  </>  </></></></> 


subroutine  mcstep 

calculate  mode  conversion  coefficients 
parameter  (maxmds-30) 
implicit  real*8  (a-h,o-z) 

common/mc  inpt/xtra(3 .maxmds) , temp(maxmds) ,freq.omega,waveno,mik, 

$  alpha , ak, akl3 ,kal3 ,ka23 ,ngsq , nthsq, skipl (3) 

common/mc  stor/a(maxmds , maxmds) , fofr (maxmds) ,s(maxmds) ,tp(maxmds) , 
$  skip2(3) .nrmode.nrslab 

common/mc  plot/iskipl(78) .nprint 

common  norm (maxmds .maxmds) , cap i (maxmds) , ans (maxmds) , ps (maxmds) , 
exO (maxmds ) ,  eyO (maxmds ) ,  hxO (maxmds ) ,  tau(maxmds) , 
pexO(maxmds) ,peyO(maxmds) ,phxO(maxmds) ,ptau(maxmds) , 
ext (maxmds ) ,  eyt (maxmds ) ,  hxt (maxmds ) ,  hyt (maxmds ) , 
pext(maxmds) , peyt(maxmds) , phxt (maxmds) ,phyt(maxmds) 
complex*16  the tap , fofr ,xtra, temp , ans , tp , a , s , ssq , sj , sm, tj , tm, ps , 
qO ,hlO ,h20 .hlprmO ,h2prm0 , qt ,hlt ,h2t .hlprmt ,h2prmt , 
al , a2 , a3 , a4 , mik , eyO j , eyOm , 

fprp , df prp , f prl , dfprl , mult , argt , argO , argg , norm , capi , 
exO,  ext,  eyO,  eyt,  hxO,  hxt,  hyt,  ngsq,  tau, 
pexO , pext , peyO , pey t , phxO , phxt , phy t , pngsq , p tau , 
one, zero, w 

real*8  kal3 ,ka23 .nthsq 
integer  pnmode 

data  one/(l . dO ,0 . dO)/, zero/(0 . dO , 0 . dO)/ 

data  w/(0.d0, -1.45749544d0)/ 

do  90  m-l.nrmode 
thetap=m 
ssq^s (m)**2 
tm-zsqrt( ngsq -ssq) 
tau(m)=tm 

a2-dcmp lx ( 0 . dO , kal 3 ) *  tm 

al=a2/ngsq 

q0-ka2  3*(one-ssq) 

call  mdhnkl(q0 ,hl0 ,h20 , hlprmO ,h2prm0 , thetap , 'MCO  ') 

exO(m)-- tm/ngsq 

ey0(m)=fofr(m) 

hxO(m)=tm*fofr (m) 

hyO  -one 

qt=ka23*(nthsq-ssq) 

call  mdhnkl(qt, hit, h2t, hlprmt, h2prmt, thetap, 'MCT  ') 

a3=hlt  *h20  -hlO  *h2t 

a4=hlt  *h2prm0-hlprm0*h2t 

fprl-a4-al*a3 

fprp-a4-a2*a3 

a3-hlprmt*h20  -hlO  *h2prmt 

a4-hlprmt*h2prm0-hlprm0*h2prmt 
dfprl=a4-al*a3 
dfprp-a4-a2*a3 

ext(m)-dcmplx(0 . dO , akl3)*dfprl/w 
eyt(ra)“eyO(m)*fprp/w 

hxt(m)-dcmplx(0 . dO , -akl3)*ey0(m)*dfprp/w 


hyt(m)“fprl/w 
90  continue 
c 

if(nrslab  .eq.  1)  then 
c 

c  first  slab 

do  32  in=l,nrmode 
do  32  j-l,nrmode 
if(j  .eq.  m)  then 
a(j  ,ni)-one 
else 

a(j ,m)=2ero 
end  if 

32  continue 

c 

else 

c 

c  integrals  in  current  slab 

if(nprint  .gt.  3)  print  906,nrslab 

do  lAO  m— 1 , nnnode 

sm=s(m) 

tm=tau(m) 

ey0m=ey0(ni) 

if(nprint  .gt.  3)  print  902 
do  140  j -1 , nrmode 
if(j  .eq.  m)  then 
c  modes  are  equal 

mult-sm* ( 2 . dO/alpha) 
ssq-sm**2 
qO-one-ssq 
qt=nthsq-ssq 

ar gt-temp (m) * ( q t*ey t (m) **2 -hxt (m) **2 ) + ( qt*hy t (ra) **2 - ext (m) **2 ) 
argO- temp (m) * ( q0*ey0 (m) **2 -hxO (m) **2 ) + (qO  - exO (m) **2 ) 

ar gg= ( one+ey 0m**2 ) *sm/ (mik*tm) 
else 

sj=s(j) 

tj-tau(j) 

eyOj-eyO(j) 

mult-one/ (mik*(sj -sm) ) 

argt=temp(m)*(eyt(m)*hxt(j)-eyt(j)*hxt(m)) 

$  -  (hyt(m)*ext(j)-hyt(j)*ext(m)) 

arg0=temp (m) * ( eyO (m) *hx0 ( j ) - eyO ( j ) *hx0 (m) ) 

^  -  (  exO(j)-  exO(m)) 

ar gg- ( one+eyO j  *eyOm) * ( s j  +sm) / (mik* ( t j  +tm) ) 
end  if 

norm(j ,m)-rault*(argt-argO) -argg 
if(nprint  .gt.  3)  print  908 ,m, j ,norm(j ,m) 

140  continue 
c 

if(nthsq  .ne.  pnthsq)  then 

c  previous  slab  had  different  topht,  must  recompute  fields 

do  200  j— l.pnmode 
the tap- j 
ssq-ps( j )**2 

a2-dcmplx(0.d0,kal3)*ptau(j) 
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al=a2/pngsq 
qO“ka2 3* ( one  - s  sq ) 

call  mdhnkl(qO ,hlO ,h20 .hlprmO ,h2prm0 , thetap , 'MCPO ' ) 
qt“ka23*(nthsq-ssq) 

call  mdhnkl(qt,hlt,h2t,hlprmt,h2prmt, thetap, 'MCPT' ) 

a3=hlt  *h20  -hlO  *h2t 

a4=hlt  *h2prmO-hlprmO*h2t 

fprl=a4-al*a3 

fprp=a4-a2*a3 

a3=hlprmt*h20  -hlO  *h2prmt 

a4=hlprmt*h2prm0  -hlprinO*h2prmt 
dfprl=a4 - al*a3 
df prp=a4 - a2*a3 

pext( j )=dcmplx(0 . dO , akl3)*dfprl/w 
peyt(j)=peyO(j )*fprp/w 

phxt( j )=dcmplx(0 . dO , -akl3)*peyO( j )*dfprp/w 
phyt(j)=fprl/w 
200  continue 

end  if 
c 

c  integrals  across  slab  boundary 

init=0 

if(nprint  .gt.  2)  print  900,nrslab 
c 

do  500  j=l,pnmode 

sj-ps(j) 

sjr=sj 

sji-=sj*(0.d0,  -l.dO) 
tj»ptau(j) 
eyOj=peyO(j ) 

if(nprint  .gt.  2)  print  902 
c 

do  400  m=l , nrmode 

sm-s(m) 

smr-sm 

snii=sin*(0 .  dO ,  - 1 .  dO) 

tffl=tau(ra) 

eyOm=eyO(ni) 

if(smr  .eq.  sjr  .and.  smi  .eq.  sji)  then 
c  modes  are  equal 

mult=sm* ( 2 . dO/alpha) 
ssq'=sm**2 
qO-one-ssq 
qt-nthsq- ssq 

argt=temp(m)*(qt*eyt(ra)**2-hxt(m)**2)+(qt*hyt(m)**2-ext(m)**2) 
argO=-temp(m)*(qO*eyO(m)**2-hxO(m)**2)  +  (qO  -ex0(m)'**2) 

argg- (one+ey0m**2 ) *sra/ (mik* tm ) 
else 

mult--one/(mik*(sj  -sra)  ) 

argt=temp(m)*(eyt(m)*phxt(j ) -peyt(j )*hxt(m) ) 

$  -  (hyt(m)*pext(j)-phyt(j)*ext(m)) 

argO-temp(ra)*(eyO(m)*phxO(j ) -peyO(j )*hxO(m) ) 

$  -  (  pexO(j)-  exO(m)) 

argg- ( one+eyO j  *eyOra) * (s j  +sm) /(mik* ( t j  +tm) ) 
end  if 
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capi(m)=-mult*(argt-argO)  -argg 
if(nprint  .gt.  3)  print  910,m,  j  ,capi(ni) 

400  continue 

c 

c  calculate  conversion  coefficients  for  hy 

call  clineq(norm,capi ,ans .nrmode .maxmds , init , err) 
init-1 

if(nprint  .gt.  3)  print  902 
do  430  m=l, nrmode 
a(j ,m)=ans(m) 
if(nprint  .gt.  2)  then 
ar-ansCm) 

ai=ans (m)*(0. dO , -1 . dO) 
db-10 . *dloglO(ar*ar+ai*ai) 
ang=datan2 (ai , ar)*57 . 2957795d0 
print  903 ,m, j ,ar ,ai,db,ang 
end  if 

430  continue 

c 

500  continue 

end  if 
c 

c  save  data  for  next  slab 

pngsq-ngsq 
pnthsq=nthsq 
pnmode-nrmode 
do  60  j-l,pnmode 
ps(j)-s(j) 
ptau(j)=tau(j) 
pex0(j)-ex0(j) 
pey0(j)-ey0(j) 
phx0(j )-hx0(j ) 
pext(j)=ext(j) 
peyt(j)=eyt(j) 
phxt(j)-hxt(j) 

60  phyt(j )-hyt(j) 
return 
c 

900  formate 'Oconversion  coefficients  for  slab  ',i2) 

901  formate  AC  ,  i2  ,  '  . '  ,  i2  ,  ’  ,  Ip2el5 . 5 , 5x, 

$  '20*logl0(A)-' ,0p2fl0.3) 

902  format('  ') 

903  formate'  le i2 i2 Ip2el5 . 5 , 5x , 

$  '20*logl0eT)-' ,0p2fl0.3) 

906  formate ' OIntegrals  in  slab',i3) 

908  formate'  Norme ' ,  i2,',',i2,')  -',lp2el5.6) 

910  formate'  Capie',  i2,'.',i2,')  -',lp2el5.6) 

end 
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subroutine  rpower 

c  subroutine  to  calculate  the  radiated  power  -  pw 

c 

parameter  (maxsgs=50) 
implicit  real*8  (a-h,o-z) 
c 

real*8  jl2,j32,j52 

real*4  pin 

real*4  gamma, phi 

real*8  iO , idl , idlsav, moment , ms 

dimension  moment (maxsgs) ,ms(maxsgs) , sbar(maxsgs) , r(maxsgs .maxsgs) 
c 

common/com01/antlen(maxsgs) , antrad(maxsgs) , antalt (maxsgs) , 

$  antang(maxsgs) ,npts 

common/com02/cpl ,psi , iclock, chicp 
common/com03/x(maxsgs) ,y(maxsgs) ,z(maxsgs) 
common/com04/xmid(maxsgs) ,ymid(maxsgs) , zmid(maxsgs) 
common/comO  5 /gamma ( maxsgs ) , phi ( maxsgs ) 
common/com06/idl (maxsgs) 

common/mc  inpt/skipl(242) ,waveno,skip2(13) 
common/mc  plot/skip3(27) ,pin, iskip4(24) 
c 

pi-3.1415926536 

nsegs=npts-l 

c 

antlen(npts+l)=antlen(npts)+cpl 
antalt (npts+l)=antalt(npts) 
c 

c  convert  lengths  to  kilometers 
do  i-l,npts+l 

antlen(i)“antlen(i)*0. 3048d-3 
antal t ( i ) -antal t ( i ) *0 . 3048d- 3 
enddo 

do  i-l,nsegs+l 

xmid(i)=xmid(i)*0 . 3048d-3 
ymid(i)=ymid(i)*0. 3048d-3 
zmid( i)— zmid( i)*0 . 3048d-3 
enddo 
c 

do  i-l,nsegs+l 

sbar(i)-antlen(i+l) -antlen(i) 

moment ( i )— sin (waveno*( (antlen(i)+antlen( i+1) )/2 . 0) ) 
ms( i)-moment( i)*sbar(i) 
enddo 
c 

sumeq-0 . 0 
sum-0 . 0 

do  n-l,nsegs+l 
do  m— 1 , n 

zp-waveno* ( zmid ( n) +zmid (m) ) 
zpsq-zp*zp 

cosgnm-cos (gamma(n) )*cos (gamma (m) ) 
ssqgnm-sin(gamma(n) )*s in (gamma (m) ) 
cospnm-cos(phi(n) -phi(ra)) 

singnm-0 . 5*sin(gamma(n) )*sin(gamma(m) )*cospnm 
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if(n  .eq.  m)  then 
c2p-cos(zp) 
szp-sin(zp) 
zpcu-zpsq*zp 

term»=(l .  0+3 . 0/zpcu*(szp-zp*czp)  )*cosgnin 
$  +(1 • 0+1 . 5/zpcu*( (1 . 0-2psq)*szp-zp*czp) )*ssqgnin 

suineq=sumeq+ms  (n)*ms  (m)*term 
else 

zm-waveno*(zniid(n)  -zniid(m)  ) 
zmsq=zm*zm 

r(n,m)-waveno*sqrt(  (xniid(n)  -xmid(m)  )**2+(ymid(n)  -yniid(m) )**2) 
rsq-r(n,m)*r(n,m) 
wm=sqrt(zmsq+rsq) 
wnmil=l .  O/wm 
winnil2=sqrt  (wmml ) 
wnim3  2“wmml*wiiunl  2 
wmmS  2“wmml*winm3  2 
wp“sqr t ( zpsq+rsq ) 
call  jfunct(win, jl2, j32,j52) 
terml-(  cosgnm+singnin)*wmml2*jl2 
term3-(  -cosgnm+singnm)*wn!m32*(  j  32-zmsq*wnmil*j  52) 
if(xmid(n)  .eq.  xmid(m))  then 
if(ymid(n)  .eq.  ymid(m))  then 
phinm-0 . 0 
else 

if(ymid(n)  .gt,  ymid(m))  then 
phinm-pi/2 . 0 
else 

phinin=-pi/2 . 0 
endif 
endif 
else 

phinm=atan2  (yinid(n)  -yinid(m)  ,xinid(n)  -xniid(m) ) 
endif 

cospm-cos(phinm-phi(ni) ) 

cospn-cos (phinm-phi(n) ) 

cos2p-cos(2 . 0*phinin-phi(n)  -phi(m)  ) 

c  s  cnni=c  o  s  ( gamma  (  n  ))*  s  in  (  gamma  (  m)  )  *c  o  spm 

cscmn-cos (gamma(m) )*sin(gamma(n) )*cospn 

ssc-0. 5*s in ( gamma (n) )*sin(gamma(m) )*cos2p 

term5-(cscnm*zm+cscmn*zm+ssc*r(n,m) )*r (n,m)*wmm52*j  52 

call  jfunct(vrp, jl2,j32,j52) 

wpml-1 . 0/wp 

wpml2-»sqrt(wpml) 

wpm3 2=wpral*wpml 2 

wpm5  2-wpml*wpm3  2 

term2— (  cosgnm-singnm)*wpml2*jl2 

term4-(  cosgnm+singnm)*wpm32*(j 32-zpsq*wpml*j 52) 

term6-( -cscnm*zp+cscmn*zp-ssc*r(n,m) )*r(n,m)*wpm52*j  52 

s\im-sum+ms (n) *ms (m) * ( terral+term2+term3 - term4+term5+term6 ) 

endif 

enddo 

enddo 

pw-20 . O*waveno*waveno*sumeq 
$  +60 . O*waveno*waveno*sum 


A-29 


iO-sqrt (pin*l . Oe3/pw) 
do  i-l,nsegs+l 

idl(i)-iO*ms(i)*l .0e3 
enddo 
return 
end 
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SUBROUTINE  JFUNCT(Q. J12 , J32 , J52) 

implicit  real*8  (a-h,o-z) 

REAL*8  J12,J32,J52 

SINQ=SIN(Q) 

COSQ-COS(Q) 

RTQ-SQRT(Q) 

J12-SINQ/RTQ 

J32-(-COSQ+Jl2/RTQ)/RTQ 

J52-3.*J32/Q-J12 

RETURN 

END 
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subroutine  tacanio 


c 

implicit  real*8  (a-h,o-z) 
parameter  (maxsgs-50 ,maxmds“30) 

c  TOWPLANE  STEADY  STATE  CONFIGURATION ... towp lane  in  circular  orbit 

c 

C  V 

c  rpl 

c  zpl 

c  amg 

c  abase 

c  cddrog 

c  amug 

c  d 

c  al 

c  picf 

c  cd 

c  X 

c  xcp 

c  alpha 

c  c appal 

c 
c 
c 

dimension  dsstor(17) 

common/com01/antlen(maxsgs) , antrad(maxsgs) , antalt(maxsgs) , 

$  antang(maxsgs) ,npts 

common/com02/cpl ,psi , iclock, chicp 
common/com09/v, rpl , zpl 

common/rac  inpt/skipl(8 .maxmds) , freq , skip2 (15) 
c 

namelist/tacin/v.rpl ,zpl,rzd,zzd,amg,psi , iclock, itrset , cpl , chicp 
data  g/32 . 17/ , pi/3 . 1416/, acl/ . 53/, acra/- . 68/ ,x/l . 34/ ,xcp/2 . 31/ , 

$  abase/5 . 6/,cddrog/. 6/,amug/l . 095e-l/,d/l . 75e-2/, dds/1 . e3/, 

$  picf/2 . 2e - 2/ , cd/1 . 03/ , delr/5 . 0/ , delz/50 . 0/ 

data  dsstor/50 . , 130 . ,4*200. ,4*300 . , 4*400 . , 2*600 . , 1000 .  / 
c 
c 

read  tacin 

print  tacin 

al“4 . 9179d+5/freq-cpl 

iter-1 

print  1101 ,v, rpl , zpl , rzd, zzd 
print  1102 , amg, abase , cddrog 
print  1103 , amug, d,al , dds , picf , cd 
ddrog-sqrt (4 . *abase/pi) 
if(cd  .le.  .1)  stop 
k-1 

istop-0 
index-0 
5  iter-iter+1 

if(iter  .gt.  30)  then 
print  1000 
stop 
endif 


towplane  true  airspeed, knots 

towplane  turn  radius{feet)  for  straight  and  level  flight 

towplane  density  altitude  ft. 

drogue  weight, lbs. 

drogue  base  area.sq.  ft. 

drogue  drag  coefficient 

towline  weight, ]b/ft 

towline  diamtcer.ft. 

towline  length,ft, 

towline  skin  friction  coefficient 

towline  drag  coefficient 

distance  from  towpoint  to  drogue  center  of  gravity, ft. 
distance  from  towpoint  to  ddrogue  center  of  pressure, ft. 
pitch  angle  of  drogue  centerline(radians) , positive  is  nose  up 
side  slip  angle  in  radians , towpoint  left  positive 
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16 


continue 

zz=zzd 

rz-rzd 

i=l 

ds“20 . 

omega-v*l . 69/rpl 
m-1 

arg-1 . - . 006875*zz/1000 . 
if(arg  .gt.  0.)then 

rhoz=0 . 002378*arg**4 .256 
else 

print  *,  'does  not  converge' 
stop 
end  if 

omegsq=omega*omega 
s=0 . 
thz-0 . 

qs-. 5*rhoz*rz*rz*omegsq*abase 
a=0 . 5*amg*x/qs+acl*x* . 04 
b- ( ac 1+cddr og) *x - ddr og*acm 
C--0 . 04*acl*x-amg*x/qs 
alpha= ( -b+sqr t (b*b - 4 . *a*c ) ) * . 5/a 
if (alpha  .ge.  .8)  then 
a2-. 85*qs*xcp 
b2“-amg*x- . 85*pi*qs*xcp 
c2“amg*x*pi*. 5- .85*qs*xcp 
alpha"( -b2 - sqrt (b2*b2 -4 . *a2*c2) ) * . 5/a2 
in-2 

am-amg/g 
end  if 

al- - amg/g*r z*omegsq - qs*acl* . 04 
bl-qs*(ddrog*acm/(x*cos (alpha) ) -acl-cddrog) 
cappal-al/bl 

if(cappal  . le .  .04)  then 
trpz--amg/g*rz*oraegsq 
else 

trpz--amg/g*rz*oinegsq+acl*(cappal- . 04)*qs 
endif 

if (alpha  .le.  .04)  then 
tzpz=amg 

trthpz=0 . 5*rhoz*abase*cddrog*rz*rz*omegsq 
else 

if(m  .eq.  1)  then 

tzpz=amg-acl*(alpha- .04)*qs 
trthpz“0 . 5*rhoz*abase*cddrog*rz*rz*omegsq 
else 

if(m  .ne.  2)  stop 

tzpz-amg- . 85*qs*sin(alpha)**2 .*cos (alpha) 
trthpz- . 85*qs*sin(alpha)**3 . 
endif 
endif 

tz-sqrt(trpz*trpz+trthpz*trthpz+tzpz*tzpz) 

rpz-trpz/tz 

thpz-tr thpz/ ( tz*rz ) 

zpz-tzpz/tz 
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rthpz=thpz*rz 

npts=l 

if (mod(index, 3)  .eq.  0  .or.  k  .eq.  4)  then 
print  1 

print  21 , s , rz , zz , thz , tz , rpz , zpz , rthpz 
endif 

if(k  .eq.  4)  then 
antlen(npts)=s 
antrad(npts)=rz 
antalt  (npts)==zz 
an tang ( np  ts ) =thz 
endif 
8  s=s+ds 

npts-npts+1 

if(npts  .gt.  maxsgs-2)  then 

print  *, 'number  of  segments  too  big' 
stop 
endif 
z=zz 
t=tz 
r=rz 
rp=rpz 
thp“thpz 
zp-zpz 
kk-1 

10  za-.5*(z+zz) 
ta-. 5*(t+tz) 
ra-. 5*(r+rz) 
thpa= . 5* ( thp+thpz ) 
zpa“.5*(zp+zpz) 
rpa-. 5*(rp+rpz) 
arg-1. - .006875*za/1000. 
if(arg  .gt.  0.)then 

rho“ . 002378*arg**4 . 256 
else 

print  *,  'does  not  converge' 
stop 
end  if 

thpasq-thpa*thpa 

rasq=ra*ra 

qd=-.  5*rho*d*rasq*omegsq*cd 
strth-sqrt(l . -rasq*thpasq) 

trp- ( ta*ra*thpasq - amug/g*ra*omegsq- qd*ra*rpa*thpa*s  tr th) *ds+trpz 
trthp-(-ta*rpa*thpa+qd*strth**3.+. 5*rho*d*picf*rasq*omegsq)*ds+ 

$  trthpz 

tzp-(amug'qd*ra*thpa*zpa*strth)*ds+tzpz 
tl-sqr t ( trp*trp+tr thp*trthp+tzp*tzp ) 
rpl-trp/tl 

rl-rz+(rpl+rpz)*. 5*ds 
thpl-tr thp/ ( tl*rl ) 
zpl-tzp/tl 

zi-zz+. 5*(zpl+zpz)*ds 
rthpl“thpl*rl 
thl-thz+(thpz+thpl)* . 5*ds 
if (abs((zl-z)/z)  .gt.  .001)  go  to  100 
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if (abs( (rl-r)/r)  .gt.  .001)  go  to  100 
if (abs ( ( tl- t)/t)  .gt.  .001)  go  to  100 
32  if(l  . le .  i  .and.  i  . le .  17)  then 
ds=dsstor(i) 
else 

if(i  .eq.  18)  then 

if (s-al+1 . Oe-4  .eq.  0.0)  then 
ds=dds 
else 

if (s-al+l.Oe-4  .It.  0.0)  then 
go  to  50 
else 

go  to  190 
endif 
endif 
else 

print  i  is  not  between  1  and  18' 

stop 
endif 
endif 
i-i+1 
50  zz-zl 
rz-rl 
tz=tl 
rpz=rpl 
thpz-thpl 
zpz”zpl 
thz-thl 
trpz-trp 
trthpz-trthp 
tzpz-tzp 

if(al-s  .It.  1000. )ds-al-s 
if(k  .eq.  4)  then 

print  21 , s , rl , zl , thl , tl , rpl , zpl , rthpl 
antlen(npts)=-s 
antrad(npts)-rl 
antalt(npts)“zl 
antang ( np  t s ) =thl 
endif 
go  to  8 

100  if(kk  .gt.  20)  then 
print  3 
go  to  32 
endif 
z-zl 
r-rl 
t=tl 
rp-rpl 
thp“thpl 
zp-zpl 
rthp-rthpl 
kk-kk+1 
go  to  10 
190  continue 

if (mod( index , 3)  .eq.  0  .or.  k  .eq.  4) 
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$  print  21 , s , rl , zl , thl , tl , rpl , zpl , rchpl 
if(k  .eq.  4)  then 
antlen(npts)=s 
antrad(npts)=rl 
antalt(npts)=zl 
antang(npts)=thl 
endif 

if^iter  .gt.  itrset)  go  to  18 
go  to(200,203,18,999) ,k 

200  if(abs((zl-zpl)/zpl)- .010)202,210,210 

202  go  to (204 , 240) , k 

203  if(abs(zl-zpl)-2000. )204,204,210 

204  if(abs((rl-rpl)/rpl)- .003)206,220,220 

206  go  to(240,200) ,k 

210  k=l 

2zd=zzd+. 6*(2pl-zl) 
go  to  5 

220  if(k  .eq.  1)  then 

drzd=. 16*(rpl-rl) 
else 

if(k  .eq.  2)  then 

if (abs( (rl'rstore)/rpl)  .le.  .002)  drzd=. 2*(rpl-rl) 
drzd=(rpl-rl)*drzd*. 6/(rl-rstore) 
if(abs(drzd)  .gt.  150.)  then 
drzd=150 .*drzd/abs (drzd) 
else 

if(abs(drzd)  .le.  5.)  drzd-'S  .  *drzd/abs  (drzd) 
endif 
else 

print  k  is  not  1  or  2’ 
stop 
endif 
endif 

rzd=rzd+drzd 

k=2 

rstore=rl 
go  to  5 

240  h=3 

18  continue 

if(index  . gt .  30)go  to  1199 
if(istop  .eq.  500)go  to  999 
if (mod(index, 3)  ,eq.  0)then 

if (abs((zl-zpl)/zpl)  .It.  .001  .and.  abs( (rl-rpl)/rpl)  .It. 

$  .001)go  to  17 

fl=zl-zpl 
f2=rl-rpl 
zlsav=zl 
rlsav-rl 
zzdsav=zzd 
rzdsav-rzd 
zzd-zzd+delz 
index- index+1 
go  to  16 
else 

if (mod(index , 3)  .eq.  l)then 
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dzldzd=(zl-zlsav)/delz 
drldzd=(rl-rlsav) /delz 
zzd=zzdsav 
rzd=rzdsav+delr 
index=index+l 
go  to  16 
else 

dzldrd=(zl-zlsav)/delr 
drldrd=(rl -rlsav)/delr 
den=dzldzd*drldrd-drldzd*dzldrd 
delzd= ( f 2*dz Idrd- f l*drldrd) /den 
delrd=(fl*drldzd-f2*dzldzd) /den 
z  z  d=z  z  ds  av+de Iz  d 
rzd=rzdsav+delrd 
index=index+l 
go  to  16 
end  if 
end  if 
17  k=4 

istop=500 
go  to  16 
999  return 

1199  print  *,'does  not  converge' 
stop 
c 

1  f ormat ( Ih  , 5x , Ihs , llx , Ihr , llx , Ihz , llx , 2hth , lOx , Iht , lOx , 2hrp , lOx , 

$  2hzp , 9x , 4hrthp) 

2  format(8el2.4) 

21  format(8el2 .4) 

3  forniat(13h  not  converge) 

1000  format ( Ihl , 14hrlnot  converge) 

1101  format (Ih  , 8hv(knot)= , el2 . 4 , 2x,4hrpl=, el2 . 4 , 2x , 4hzpl= , el2 . 4 , 2x, 

$  4hrzd= , el2 . 4 , 2x,4hzzd=, el2 . 4) 

1102  format (Ih  , 4hamg=- , el2 . 4 , 2x, 6habase=, el2 . 4 , 2x , 7hcddrog= , el2 . 4) 

1103  format(lh  , 5hamug= , el2 .4 , 2x, 2hd“, el2 . 4 , 2x, 3hal=, el2 . 4 , 2x. 4hdds=, 

$  el2 . 4 , 2x , 5hpicf=,el2 .4 , 2x, 3hcd= , el2 . 4) 

1203  format(lh  , 3x, 6halpha= , el2 .4, 5x, 7hcappal=, el2 . 4) 
c 

end 
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SUBROUTINE  HTGAIN ( lOPT . FREQ , SIGMA , EPSR , ALPHA , NRMODE , TP . Z , HG ) 
IMPLIGIT  REALMS  (A-H,0-Z) 

COMPLEX*16  TP(1) ,HG(3,1) .HGO, 

$  C , S  SQ , NGSQ , SQROOT . RATIO , A1 , A2 , A3 , A4 , MI , ONE 

$  PO , HIO , H20 , HlPRMO , H2PRM0 , HIZ , H2Z , HIPRMZ , H2PRMZ , EXPZ 

REALMS  K,KA13,KA23 

DATA  MI/(0.D0, -1. DO)/, ONE/(l. DO, 0. DO)/, HG0/(0. DO, 1 .4574954400)/ 
DATA  DTR/1.745329252D-02/ 

C 

NGSQ=DCMPLX(EPSR, -SIGMA/(5.5633459D-8*FREQ)) 

K=2.0958426D-02*FREQ 

IF(ALPHA  .EQ.  O.DO)  GO  TO  5 

AK=ALPHA/K 

AK13=DEXP(DL0G(AK) /3 .DO) 

AK23=AK13**2 
KA13=1.D0/AK13 
KA23=KA13**2 
P1=KA23*ALPHA*Z 
EXPZ=DEXP( . 5D0*ALPHA*Z) 

5  DO  20  M=l, NRMODE 

SSQ=ZSIN(TP(M)*DTR)**2 

CSQ=ONE-SSQ 

SQROOT=ZSQRT(NGSQ-SSQ) 

IF(DIMAG(TP(M))  .LE.  -10. DO  .OR.  ALPHA  .EQ.  O.DO)  GO  TO  10 
PO=KA23*(ONE-SSQ) 

CALL  MDHNKL(PO  , HIO, H20, HlPRMO, H2PRM0, TP(M) , 'HG  1') 

GALL  MDHNKL(P0+P1, HIZ, H2Z, HIPRMZ, H2PRMZ, TP(M) , 'HG  2') 

Al-HlO  *H2Z  -HIZ  *H20 

A2=H1PRM0*H2Z  -HIZ  *H2PRM0 

A3=H10  *H2PRMZ-H1PRMZ*H20 
A4=H1PRM0*H2  PRMZ - H 1 PRMZ*H2  PRMO 
RATIO=SQROOT/NGSQ 
G- . 5DO*AK23+KA13*MI*RATIO 
HG ( 1 , M) =EXPZ* ( C*A1+A2 ) 

HG  (  2  ,  M)  =KA1 3*MI*SQROOT*Al-i-A2 

HG ( 3 , M) - . 5D0*AK*MI*HG ( 1 , M) +AK13*MI*EXPZ* ( C*A3+A4 ) 

IF(IOPT  .EQ.  1)  GO  TO  20 
HG(1,M)=HG(1,M)/HG0 
HG(2,M)=HG(2,M)/HGO 
HG(3,M)=HG(3,M)/(RATIO*HGO) 

GO  TO  20 

10  C-ZSQRT(ONE-SSQ) 

EXPZ=ZEXP(DCMPLX(0 . DO ,K*Z)*C) 

Al-=  (  NG  SQ*C  -  SQROOT )  /  ( NGSQ*G+SQROOT  ) 

A2=(C- SQROOT)/ (G+SQROOT) 

HG ( 1 , M) =EXPZ+Al/EXPZ 
HG(2 ,M)-EXPZ+A2/EXPZ 
HG(3 ,M)=(EXPZ-A1/EXPZ)*G 
IFCIOPT  .EQ.  1)  GO  TO  20 
HG(1,M)-HG(1,M)/(0NE+A1) 

HG(2,M)-HG(2,M)/(ONE+A2) 

HG(3,M)-HG(3,M)/((0NE-A1)*C) 

20  CONTINUE 
RETURN 
END 


A-38 


-o>-c/>-c/>-a>-c/></> -»/></>-(/></></></>■(/></>  </></><j></></><j><rK/></></></><n</><j>  -c/></></>-c/>-</>-c/>-c/>-<y>-c/>-c/>-c/></>-c/>-</> 


SUBROUTINE  MDHNKL  (Z , HI , H2 .HIPRME , H2PRME , THETA , IDBG) 

IMPLICIT  COMPLEX*16  (A-H,0-2) 

COMPLEX*! 6  I , MPOWER , MTERM 

REAL*8  A,B,C,D,CAP,PART1,PART2,ZMAG 

CHARACTER*^  IDBG 

DIMENSION  A(30),  B(30) ,  C(30),  D(30) ,  CAP(30) ,  PART1(2),  PART2(2) 
EQUIVALENCE  (PARTI , TERM4) ,  (PART2,SUM4) 

DATA  A  /  9.3043671692922944819D-01,  3 . 1014557230974314911D+01 , 

2 . 0676371487316209897D+02 ,  5 . 7434365242545027449D+02 , 

8. 7021765519007617234D+02,  8 . 2877871922864397320D+02 , 

5 . 4168543740434246542D+02 ,  2 . 5794544638302022111D+02 , 

9.3458495066311674231D+01,  2 . 6626351870744066662D+0I , 

6.1210004300561072794D+00,  1 . 1592803844803233472D+00 , 

1 . 8401275944132116616D-01 ,  2 . 4833030963741048003D-02 , 

2.8842080097260218300D-03,  2 . 9133414239656786138D-04 , 

2.5827494893312753646D-05,  2 . 0256858739853140063D-06 , 

1 . 4155736366074870734D-07 ,  8 . 8695090013000443124D-09 , 

5.0110220346327933889D-10,  2 . 5658074934115685526D-11 , 

1.1961806496091228666D-12.  5 . 0988092481207283185D- 14 , 

1.9948392989517716388D-15,  7 . 1886100863126905797D- 17 , 

2.3938095525516785112D-18,  7 . 3883010881224645255D-20, 

2.1194208514407528762D-21,  5 . 6653858632471341093D-23/ 

DATA  B  /  6.7829872514427588456D-01,  1 . 1304978752404598033D+01 , 

5 . 3833232154307609704D+01 ,  1 . 1962940478735024376D+02 , 

1 . 5337103177865415841D+02 ,  1 . 2780919314887846509D+02 , 

7.4742218215718400631D+01,  3 . 2355938621523117060D+01 , 

1 . 0785312873841039006D+01 ,  2 . 8532573740320209005D+00 , 

6.1360373635097223595D-01,  1 . 0937678009821251966D-01 , 

1.6422939954686564465D-02,  2 . 1055051223957133911D-03 , 

2.3316778764072130571D-04,  2 . 2528288660939256561D-05 , 

1.9156708045016374595D-06,  1.4446989475879618839D-07, 

9.7286124416697769730D-09,  5 . 8854279743918795891D- 10 , 

3 . 2160808603234314644D-11 ,  1 . 5952782045255116351D-12 , 

7.2151886229105003778D-14,  2 . 9876557444763976717D- 15 , 

1.1368553061173507104D-16,  3 . 9889659863766691603D-18 , 

1.2946984700995355913D-19,  3 . 8985199340546088228D- 21 , 

1.0920223904914870636D-22,  2 . 8527230681595795812D-24/ 

DATA  C  /  4.6521835846461472410D-01,  6 . 2029114461948629822D+00 , 

2 . 5845464359145262382D+01 ,  5 . 2213059311404570392D+01 , 

6.2158403942148298012D+01,  4 . 8751689366390821897D+01 , 

2.7084271870217123228D+01,  1 . 1215019407957400909D+01 , 

3.5945575025504490022D+00,  9 . 1815006450841609147D-01 , 

1.9128126343925335199D-01,  3 . 3122296699437809740D-02 , 

4.8424410379295043444D-03,  6 . 0568368204246458321D-04 , 

6.5550182039227768583D-05,  6 . 1985987743950608612D-06 , 

5 . 1654989786625507119D-07 ,  3 . 8220488188402150986D-08 , 

2.5278100653705126277D-09,  1 . 5033066103898380141D-10 , 

8.0822936042464409157D-12,  3 . 9473961437101054471D- 13 , 

1.7590891906016512675D-14,  7 . 1814214762263778920D-16 , 

2.6957287823672589641D-17,  9 . 3358572549515461865D- 19 , 

2.9922619406895981315D-20,  8 . 9015675760511620701D- 22 , 

2.4644428505125033375D-23,  6 . 3656020935361057409D-25/ 

DATA  D  /  6.7829872514427588456D-01,  4 . 5219915009618392131D+01 , 

$  3.7683262508015326776D+02,  1 . 1962940478735024344D+03 , 

$  1.9938234131225040548D+03,  2 . 0449470903820554375D+03 , 
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1.4201021460986496090D+03,  7 . 1183064967350857463D+02 , 

2 . 6963282184602597492D+02 ,  7 . 9891206472896585111D+01 , 

1.9021715826880139294D+01,  3 . 7188105233392256682D+00 , 

6.0764877832340288572D-01,  8 . 4220204895828535644D-02 , 

1.0026214868551016149D-02,  1 . 0363012784032058021D-03 , 

9 . 3867869420580235442D-05 ,  7 . 5124345274574017960D-06 , 

5.3507368429183773360D-07,  3 . 4135482251472901638D-08 . 

1.9618093247972931935D-09,  1 . 0209780508963274472D- 10 , 

4.8341763773500352579D-12,  2 . 0913590211334783723D-13 , 

8.2990437346566602039D-15,  3 . 0316141496462685641D- 16 , 

1.0228117913786331176D-17,  3 . 1967863459247792364D-19 , 

9 . 2821903191776400453D-21,  2 . 51039629998043003090-22/ 

CAP  /  1. 04166666666666666630-01,  8.35503472222222221160-02, 

1 . 28226574556327160190-01 ,  2 . 91849026464140463150-01 , 

8.81627267443757648740-01,  3.32140828186276752640+00, 

1.49957629868625545460+01,  7.89230130115865175300+01, 

4 . 74451538868264318870+02 ,  3 . 20749009089066190040+03 , 

2.40865496408740046050+04 ,  1 . 98923119169509791210+05 , 

1 . 79190200777534380630+06 ,  1 . 74843771800341210230+07 , 

1 . 83707379676330729780+08 ,  2 . 06790403294515515080+09 , 

2.48275193759358884720+10,  3.16694549817348873150+11, 

4.27711268651347155820+12,  6.09711324113925607490+13, 

9.14869422343563967920+14,  1.44135251700093501010+16, 

2 . 37888443951757579420+17 ,  4 . 10460816009469218850+18 , 

7 . 39000494157048539930+19 ,  1 . 38592200046039431410+21 , 

2 . 70308259302757616230+22 ,  5 . 47474786196455733350+23 , 

1 . 14989370143863335240+25 ,  2 . 50141806927536039690+26/ 

OATA  1/(0.00,1.00)/ 

OATA  ONE/(l.O0,0.O0)/,TWO/(2.D0,0.D0)/,ZERO/(0.O0,0.D0)/ 

OATA  R00T3/(1. 7320508075688800,0.00)/ 

OATA  ALPHA/(8. 536672188389510-1,0.00)/ 

OATA  CONSTl/(  2.588190451025220-01,-9.659258262890670-01)/ 

DATA  C0NST2/(  2.588190451025220-01,  9.659258262890670-01)/ 

DATA  CONST3/(-9. 659258262890670-01,  2.588190451025220-01)/ 

DATA  C0NST4/(- 9. 659258262890670-01, -2. 588190451025220-01)/ 

C 

ZPOWER-ONE 

SUM3=ZER0 

SUM4=ZER0 

ZMAG=ZABS(Z) 

IF(ZMAG  .GT.  6.100)  GO  TO  70 

SUM1=ZER0 

SUM2»ZER0 

ZTERM--Z**3/(200.D0,0.D0) 

DO  50  M-1,30 

SUM1-SUM1+DCMPLX(A(M) , 0 . DO)*ZPOWER 
SUM2-SUM2+DCMPLX(B(M) , 0 . DO)*ZPOWER 
SUM3-SUM3+DCMPLX(C(M) ,0 . DO)*ZPOWER 
TERM4-DCMPLX(D(M) , 0 . DO)*ZPOWER 
SUM4-SUM4+TERM4 

IF(DABS(PART1(1)/PART2(1))  .LE.  1.0-17  .AND. 

$  DABS(PART1(2)/PART2(2))  .LE.  1.0-17)  GO  TO  60 
50  ZPOWER-ZPOWER*ZTERM 

60  GM2F-I*(Z*SUM2-TWO*SUMl)/ROOT3 

GPMFP-I* ( SUM4+TWO*Z*Z*SUM3 ) /R00T3 
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H1=Z*SUM2+GM2F 
H2=H1-TW0*GM2F 
H1PRME=SUM4+GPMFP 
H2PRME=H1PRME-TW0*GPMFP 
GO  TO  999 
70  MPOWER=ONE 

SUM1=0NE 
SUM2=ONE 

c  RTZ=CDSQRT(Z) 

RTZ=ZSQRT(Z) 

SQRTZB=RTZ*Z 

ZTERM=I/SQRTZB 

MTERM=-ZTERM 

DM=ZERO 

TERK3=ONE 

DO  80  M=l,30 

ZPOWER=ZPOWER*ZTERM 

MPOWER-MPOWER*MT  ERM 

DM-DM+ONE 

TERK1=DCMPLX(CAP(M) , 0 . DO)*ZPOWER 

TERM2=DCMPLX(CAP(M) , 0 . DO ) EMPOWER 

IF(ZABS(TERM2/TERM3)  .GE.  l.DO)  GO  TO  81 

SUMl-SUMl+TERMl 

SUM2=SUM2+TERM2 

SUM3=SUM3+DM*TERM1 

TERM4=DM*TERM2 

SUM4=SUM4+TERM4 

IF(DABS(PART1(1)/PART2(1))  .LE.  l.D-17  .AND. 

$  DABS (PARTI (2)/PART2 (2))  .LE.  l.D-17)  GO  TO  81 

80  TERM3=TERM2 

81  ZTERN-(-1.5D0,0.D0)/Z 
SUM3=SUM3*ZTERM 
SUM4-SUM4*ZTERM 

TERMl-( ( -  0 . 25D0 , 0 . DO) - I*SQRTZB)/Z 
TERM2=( ( -0 . 25D0 , 0 . DO)+I*SQRTZB)/Z 
c  EXP1=CDEXP( (0 . DO , 0 . 666666666666666667D0)*SQRTZB) 

EXP1-ZEXP( (0 . DO , 0 . 666666666666666667D0)*SQRTZB) 
EXP2=C0NST1*EXP1 
EXP3=C0NST2/EXP1 
EXP4=C0NST3*EXP1 
EXP5=C0NST4/EXP1 
c  ZTERM=ALPHA/CDSQRT(RTZ) 

ZTERM=ALPHA/ZSQRT ( RTZ ) 

TERM4-Z 

IF(PART1(1)  .GE.  O.DO  .OR.  PART1(2)  .GE.  O.DO)  GO  TO  90 
H1=ZTERM*(EXP2*SUM2+EXP5*SUM1) 

HI PRME=ZTERM* ( EXP2* ( SUM2*TERM2+SUM4 ) +EXP5* ( SUM1*TERM1+SUM3 ) ) 
GO  TO  110 

90  H1-ZTERM*EXP2*SUM2 

H1PRME-ZTERM*EXP2* ( SUM2*TERM2+SUM4 ) 

110  IF(PART1(1)  .GE.  O.DO  .OR.  PART1(2)  .LT.  O.DO)  GO  TO  120 
H2-ZTERM*  ( EXP3*SUM1+FJCP4*SUM2  ) 

H2  PRME-ZTERM* ( EXP3* ( SUM1*TERM1+SUM3 ) +EXP4* ( SUM2*TERM2+SUM4 ) ) 
GO  TO  999 

120  H2-ZTERM*EXP3*SUM1 
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H2PRME=ZTERM*EXP3* ( SUM1*TERM1+SUM3 ) 

C  CALCULATE  WRONSKIAN  AS  PARTIAL  CHECK  ON  VALIDITY 

999  SUM4=H1*H2PRME-H1PRME*H2 
IF(DABS(PART2(1))  .LE.  l.D-8  .AND. 

$  DABS(PART2(2)+1.457495441040461D0)  .LE.  l.D-8)  GO  TO  1000 

PRINT  1001, SUM4, THETA, IDBG 

1000  RETURN 

1001  FORMAT('  *****  POSSIBLE  ERROR  IN  MDHNKL:  W  =  ',1P2E15.6, 

$  '  FOR  THETA  =  ',0P2F10.4,'  AT  ' ,A4) 

END 
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SUBROUTINE  CLIN  EQ  (A,  B,  X,  N,  N  DIM,  IFLAG,  ERR) 

C 

CLIN  EQ  USES  L-U  DECOMPOSITION  TO  FIND  THE  TRIANGULAR  MATRICES  L 
AND  U  SUCH  THAT  L*U=A.  THE  MATRICES  L  AND  U  ARE  STORED  IN  A. 

THIS  FORM  IS  USED  WITH  BACK  SUBSTITUTION  TO  FIND  THE  SOLUTION  X  OF 
A*X=L*U*X=B.  N  IS  THE  NUMBER  OF  EQUATIONS  AND  NDIM  IS 
THE  DIMENSION  OF  ALL  ARRAYS.  ERR  IS  THE  ESTIMATED  RELATIVE  ERROR 
OF  THE  SOLUTION  VECTOR. 

IF  IFLAG=0,  THEN  L,  U  AND  X  ARE  COMPUTED.  OTHERWISE,  IT  IS 
ASSUMED  THAT  L  AND  U  HAVE  BEEN  COMPUTED  IN  A  PREVIOUS  CALL  AND  ARE 
C  STILL  STORED  IN  A. 

C 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEX*16  A,  B,  X,  T 

DIMENSION  A(N  DIM,  N  DIM),  B(NDIM) ,  X(N  DIM) 

DIMENSION  IROW(51),  Q(51) 

DATA  EPS  /l.D-15/ 

C 

IF  (N  .GT.  51)  THEN 
PRINT  999 
STOP 
END  IF 

IF  (N  .EQ.  1)  THEN 
X(1)=B(1)/A(1,1) 

RETURN 
END  IF 

IF  (IFLAG  .EQ.  0)  THEN 
DO  50  1=1, N 

Q(I)=0.D0 
DO  40  J=1,N 

QQ-ZABS(A(I,J)) 

40  IF  (Q(I)  .LT.  QQ)  Q(I)=QQ 
IF  (Q(I)  .EQ.  O.DO)  THEN 
LOOP=40 
GO  TO  900 
END  IF 

50  CONTINUE 
ERR=EPS 
PPIV=0 . DO 
DO  100  1=1, N 

100  IROW(I)=I 
C 

DO  500  L=1,N 

PIVOT=0.D0 

K-L-1 

DO  240  I-L,N 
IF  (K  .GE.  1)  THEN 
DO  220  J-1,K 

220  A(I,L)-A(I,L)-A(J,L)*A(I,J) 

END  IF 

F-ZABS(A(I,L))/Q(I) 

IF  (PIVOT  .LE.  F)  THEN 
PIVOT-F 
NPIVOT-I 
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END  IF 

240  CONTINUE 

IF  (PIVOT  .EQ.  O.DO)  THEN 
LOOP=240 
GO  TO  900 
END  IF 

IF  (PPIV  .GT.  PIVOT)  ERR=ERR*PPIV/PIVOT 
PPIV-PIVOT 

IF  (NPIVOT  .NE.  L)  THEN 
Q(NPIVOT)=Q(L) 

J=IROW(L) 

IROW(L)=IROW(NPIVOT) 

IROW(NPIVOT)=J 
DO  260  I-1,N 

T=A(L,I) 

A(L,I)=A(NPIVOT,I) 

A(NPIVOT,I)=T 
260  CONTINUE 
END  IF 

IF  (L  .NE.  N)  THEN 
*  T-(1.D0,0.D0)/A(L,L) 

K-L+1 

M-L-1 

DO  450  I=K,N 
IF  (M  .GE.  1)  THEN 
DO  350  J-1,M 

350  A(L,I)-A(L,I)-A(L,J)*A(J,I) 

END  IF 

A(L,I)=T*A(L,I) 

450  CONTINUE 
END  IF 

500  CONTINUE 

IF  (ERR  .GT.  l.D-5)  PRINT  998, ERR 
END  IF 

DO  620  1-2, N 

620  X(I)=(0.D0,0.D0) 

J-IROW(l) 

X(1)=B(J)/A(1,1) 

DO  700  1=2, N 

J=IROW(I) 

K-I-1 

DO  650  L=1,K 

650  X(I)-X(I)+A(I,L)*X(L) 
X(I)-(B(J)-X(I))/A(I,I) 

700  CONTINUE 
K-N-1 

DO  800  I-1,K 

J=N-I 

M-J+1 

DO  800  L-M,N 
X(J)-X(J)-X(L)*A(J,L) 

800  CONTINUE 
RETURN 

900  PRINT  997, LOOP 
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STOP 

997  FORMAT('0******ERROR  IN  CLIN  EQ,  MATRIX  IS  SINGULAR,  AT  LOOP  ',13) 

998  FORMAT ( '  0****-**CAUTION- '  , 

$  'CLIN  EQ  HAS  DECOMPOSED  AN  ILL-CONDITIONED  MATRIX. '/15X, 

$  'RESULTS  WILL  HAVE  RELATIVE  ERROR  =',1PE12.2) 

999  FORMAT('0******ERROR  IN  CLIN  EQ,  MATRIX  SIZE  GREATER  THAN  51') 

END 
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subroutine  pltbgn 
character*!  answr 
logical  first, autopl 
dimension  ia(8) , ibuff (2) 
data  ia/82,79,57,48.73,87,73,80/ 
c  ASCII  R0901WIP 

data  first/. true./, autopl/. false./ 
if(first)  then 

open(unit=52  ,  f  ile='  /<iev/plt7550a'  ) 

print  *,'If  this  is  the  hp  7550  plotter  and  you  want  auto  feed,' 
print  *,'then  set  up  the  plotter,  load  a  sheet  and  answer  y: ' 
print  *,'Do  you  want  auto  feed?' 
read  1 ,  answr 
1  format(al) 

if (answr  .eq.  'y'  .or.  answr  .eq.  'Y' )  then 
autopl=. true . 
else 

autopl= . false . 
end  if 
end  if 

if ( .not. autopl  .or.  first)  then 

print  *,'Set  up  plotter,  enter  rotation  (y/n)  when  ready' 
read  1 ,  answr 
first-. false . 
end  if 

call  hpinit(2 , 0 , 0 , 0 , 52) 

if (answr  .eq.  'y'  .or.  answr  .eq.  'Y') 

$  call  buff (1, ia.xbuff ,8) 
call  newpen(l) 
return 
end 

subroutine  plton 
dimension  ia(3) 
data  ia/27,46,89/ 
c  ASCII  esc  .  Y 

call  buff (1 , ia,xbuff , -3) 

call  newpen(l) 

return 

end 

subroutine  pltend 

call  newpen(O) 

entry  pi toff 

call  plot(0. 0,0. 0,999) 

return 

end 
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subroutine  bordr2 (xlng.xmin.xmax.xticl , xtic2 , ndx , 

$  ylng,ymin,ymax,yticl ,ytic2 ,ndy) 

c 

xscale=xlng/(xmax-xmin) 

yscale=ylng/(ymax-yTiiin) 

if (xticl*xscale  . le .  0.  .or.  xtic2*xscale  .1e.  0.)  go  to  999 
if (yticl*yscale  .le.  0.  .or.  ytic2*yscale  .le.  0.)  go  to  999 
c 

if(iabs(ndx)  .gt.  9)  then 

SX-.  15 

nx-ndx- (ndx/10)*10 
else 
sx= .  1 
nx=ndx 
end  if 
xo= . 5*sx 

if(iabs(ndy)  .gt.  9)  then 
sy= . 15 

ny-ndy- (ndy/10)*10 
else 
sy-.  1 
ny=ndy 
end  if 
yo-. 5*sy 
c 

xres”=abs  (xticl)/10 . 
yres-abs(yticl)/10. 
c 

tl-.05 

t2-.10 

yval=yinin 

ytc2“yinin 

xp”0 . 

yp=0. 

go  to  115 

112  yval—yval+yticl 

if (abs(yval-yinax)  .le.  yres)  then 
call  plot(0 . ,ylng, 2) 
if (abs(yval-ytc2)  .le.  yres)  then 
xln=- sy*(3+ny) 
yln=ylng-yo 

if(abs(yval)  .gt.  yres  .or.  abs(yval)  .ge.  10.) 

$  xln=xln-sy*aint(aloglO(abs(yval) ) ) 

if(abs(yval)  .It.  yres)  yval=-0. 
if(yval  .It.  0.)  xln=xln-sy 
call  plot(xln,yln, 3) 
call  number (xln, yin, sy, yval, 0. ,ny) 
end  if 

call  plot(0. ,ylng,3) 
go  to  120 
end  if 

yp-(yval-ymin)*yscale 
call  plot(xp,yp,2) 

if (abs(yval-ytc2)  .gt.  yres)  go  to  118 
call  plot(t2,yp,2) 
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xln=-sy*(3+ny) 
yln=yp-yo 

if(abs(yval)  .gt.  yres  .or.  abs(yval)  .ge.  10.) 
$  xln=xln-sy*aint(aloglO(abs(yval) ) ) 
if(abs(yval)  .It.  yres)  yval=0. 
if(yval  .It.  0.)  xln=xln-sy 
call  plot(xln,yln, 3) 
call  number (xln, yin, sy ,yval , 0 . ,ny) 
ytc2=ytc2+ytic2 
go  to  119 

118  call  plot ( tl , yp , 2) 

119  call  plot(xp,yp, 3) 

if (abs(yval-ymax)  .gt.  yres)  go  to  112 
c 

120  yp=ylng 
tl=ylng- .05 
t2=ylng- . 10 
xval-=xmin 
xtc2— xrain+xtic2 

122  xval-xval+xticl 

if (abs(xval-xmax)  .gt.  xres)  go  to  123 
call  plot(xlng,ylng, 2) 

if (abs(xval-xtc2)  .le.  xres)  xtc2=xtc2+xtic2 
go  to  130 

123  xp=(xval-xmin)*xscale 
call  plot(xp ,yp , 2) 

if (abs(xval-xtc2)  .gt.  xres)  go  to  128 
call  plot(xp , t2 , 2) 
xtc2=xtc2+xtic2 
go  to  129 

128  call  plot(xp , tl , 2) 

129  call  plot(xp ,yp , 2) 

if (abs(xval-xmax)  .gt.  xres)  go  to  122 
c 

10  xp-xlng 

tl=xlng- . 05 
t2=-xlng- .  10 
ytc2=ytc2-ytic2 

if (abs(ytc2-ymax)  . le .  yres)  go  to  135 

132  yval=yval-yticl 

if (abs (yval-ymin)  .gt.  yres)  go  to  133 
call  plot(xlng, 0. , 2) 
go  to  140 

133  yp=(yval-ymin)*yscale 
call  plot(xp,yp, 2) 

if (abs(yval-ytc2)  .gt.  yres)  go  to  138 
call  plot(t2,yp,2) 

135  ytc2-ytc2-ytic2 
go  to  139 

138  call  plot(tl,yp,2) 

139  call  plot(xp ,yp , 2) 

if (abs(yval-3niiin)  .gt.  yres)  go  to  132 
c 

140  yp”0. 
tl-.05 


A-48 


t2=.10 
yln=- 2 . *sx 
xtc2=xtc2 -xtic2 

if (abs (xtc2 -xmax)  .le.  xres)  go  to  145 
142  xval=xval -xticl 

if (abs(xval-xmin)  . le .  xres)  then 
call  plot(0 . , 0 . ,2) 
if (abs(xval-xtc2)  . le .  xres)  then 
xln=-xo*(2+nx) 

if(abs(xval)  .gt.  xres  .or.  abs(xval)  . ge .  10.) 

xln=xln-xo*aint(aloglO(abs (xval) ) ) 
if(abs(xval)  .It.  xres)  xval=0. 
if (xval  .It.  0.)  xln=xln-xo 
call  plot(xln,yln, 3) 
call  number (xln , yin, sx , xval , 0 . ,nx) 
end  if 

call  plot(0. ,0. ,3) 
return 
end  if 

xp=(xval-xmin)'>‘'xscale 
call  plot(xp,yp, 2) 

if (abs (xval -xtc2)  .gt.  xres)  go  to  148 
call  plot(xp , t2 , 2) 

145  xln=xp-xo*(2+nx) 

if(abs(xval)  .gt.  xres  .or.  abs(xval)  .ge.  10.) 

$  xln=xln-xo*aint(aloglO(abs(xval) ) ) 
if(abs(xval)  .It.  xres)  xval=0. 
if (xval  .It.  0.)  xln=xln-xo 
call  plot(xln,yln, 3) 
call  number (xln, yin, sx, xval, 0. ,nx) 
xtc2=xtc2 -xtic2 
go  to  149 

148  call  plot(xp , tl , 2) 

149  call  plot(xp ,yp , 3) 

if (abs(xval-xmin)  .gt.  xres)  go  to  142 
c 

999  print  100 , xlng ,xmin, xmax , xticl ,xtic2 ,ndx, 

$  ylng,ymin,ymax,yticl ,ytic2 ,ndy 

100  format ('OError  in  B0RDR2;'/ 

$  'Oxlng,  xmin,  xmax,  xticl,  xtic2,  ndx  ■=  ’ , lp5ell . 3 , 15/ 

$  'Oylng,  ymin,  yroax,  yticl,  ytic2,  ndy  =  ' , lp5ell . 3 , i5) 

call  pltend 
stop 
end 
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subroutine  curve (x , y , up , nrpts , xmin , ymin , xinc , y inc , 1 ine ) 
c 

c  x,y,up  must  be  dimensioned  at  least  nrpts 
c  xmin, ymin  are  x,y  origin  in  user  units 
c  xinc.yinc  are  x,y  scales  in  user  units  per  inch 
c 

c  line=l:  solid 

c  2 :  long  dash 

c  3 :  medium  dash 

c  4:  short  dash 

c  5 ;  dotted 

c  6:  short  +  long  dash 

c  7 :  short  +  short  +  long  dash 

c 

logical  up,upl,up2 

dimension  ipen(8) , joc(7) ,x(nrpts) ,y(nrpts) ,up(nrpts) 
data  ipen/2,2,2.3,2,3,2.3/.joc/18,  11,  14,  23,  32,  41,  16/ 
data  delr/.l/ 
c 

if (nrpts  .le.  1)  go  to  99 
c 

if(line)  1,2,3 

1  kk=»mod(line ,  7)+7 
go  to  4 

2  kk=0 

go  to  4 

3  kk~mod(line , 7) 

4  kk=kk+l 
jo-joc(kk)/10 
jc-joc(kk) -10*jo 
ip-ipen(jo) 

c 

j-0 
dr-O. 
rhol-0 . 
rho2-delr 

pxl=(x(l) -xrain)/xinc 
pyl“(y(l) -ymin)/yinc 
upl-up(l) 

if (.not.  upl)  then 
c 

c  go  to  first  position  with  pen  up 
call  plot(pxl ,pyl , 3) 
if(kk  .eq.  6)  then 
px2-(x(2) -xrain)/xinc 
py2-(y(2) -ymin)/yinc 
delx-px2-pxl 
dely-py2-pyl 

rho-s  qr t ( de lx**2+de ly**2 ) 
if(rho  .eq.  0.)  then 
dx  6-delx*.l 
dy  6-dely*.l 
else 

dx  6-delx*delr/rho*. 1 
dy  6-dely*delr/rho*. 1 
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end  if 

call  plot(pxl+dx6 ,pyl+dy6 , 2) 
end  if 
end  if 
c 

do  40  i=2,nrpts 
px2“(x(i) -xmin) /xinc 
py2=(y(i) -ymin)/yinc 
up2=up(i) 
if(up2)  then 
dr-0. 
rhol-0. 
rho2=delr 
go  to  39 
end  if 
if(upl)  then 

c  pen  has  been  up,  prepare  to  lower  pen 
call  plot (px2 , py2 , 3) 
go  to  39 
end  if 

if(kk  .eq.  2)  go  to  38 

delx-»px2-pxl 

dely->py2-pyl 

rho“sqrt(delx**2+dely**2) 

rhol-rhol+rho 

if(rho2  .gt.  rhol)  go  to  38 

delx»=delx*delr/rho 

dely-dely*delr/rho 

dx  6**delx*.l 

dy  6-=dely*.l 

if(dr  .eq.  0.)  go  to  20 

dx“delx*dr/delr 

dy=dely*dr/delr 

pxl=pxl+dx 

pyl=pyl+dy 

go  to  21 

20  if(rho2  .gt.  rhol)  go  to  38 

pxl-pxl+delx 
pyl-pyl+dely 

21  call  plot(pxl ,pyl , ip) 

if(kk  .eq.  6)  call  plot(pxl+dx6 , pyl+dy6 , 2) 

j-J+1 

ip-ipen(jo+mod(j ,jc)) 
rho2“rho2+delr 
go  to  20 

38  call  plot(px2 ,py2 , ip) 
dr“rho2-rhol 

39  pxl-px2 
pyl-py2 
upl— up2 

40  continue 

99  return 

end 
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APPENDIX  B—PROGRAM  LISTING  OF  “TACAMO 
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c  TACAMO:  TOWLINE  CONFIGURATION  CODE  FOR  TOWPLANE  IN  STEADY  STATE 
c  CONFIGURATION. .. towplane  in  circular  orbit. .. (HUANG, NADC- 

c  AM-6849, June  1969) 

c 


c 

V 

towplane  true  airspeed, knots 

c 

rpl 

towplane  turn  radius(feet)  for  straight  and  level  flight 

c 

zpl 

towplane  density  altitude  ft. 

c 

amg 

drogue  weight, lbs. 

c 

abase 

drogue  base  area,sq.  ft. 

c 

cddrog 

drogue  drag  coefficient 

c 

amug 

towline  weight, Ib/ft 

c 

d 

towline  diameter, ft. 

c 

al 

towline  length, ft. 

c 

picf 

towline  skin  friction  coefficient 

c 

cd 

towline  drag  coefficient 

c 

X 

distance  from  towpoint  to  drogue  center  of  gravity, ft. 

c 

xcp 

distance  from  towpoint  to  ddrogue  center  of  pressure, ft. 

- 

alpha 

pitch  angle  of  drogue  centerline(radians) , positive  is  nose  up 

c 

cappal 

side  slip  angle  in  radians , towpoint  left  positive 

c 

c 

namelis t/datum/v , rpl , zpl , al , rzd, zzd , amg , itrset 
c 

dimension  dsstor(17) 
c 

data  g/32 . 17/ , pi/3 . 1416/ , acl/ . 53/ , acm/- . 68/ , x/1 . 34/ , xcp/2 . 31/ , 

$  abase/5 . 6/,cddrog/. 6/,amug/l . 095e-l/, d/1 . 75e-2/, dds/1 . e3/, 
$  picf /2 . 2e -  2/ , cd/1 . 03/ , delr/5 . / , delz/50 . / 

data  dsstor/50. ,130. ,4*200. ,4*300. ,4*400. ,2*600. ,1000./ 
c 

read  datiom 
iter=l 

print  1101 ,v, rpl , zpl , rzd, zzd 
print  1102 , amg, abase .cddrog 
print  1103 , amug, d, al , dds , picf , cd 
ddrog-sqrt(4 . *abase/pi) 
if(cd  .le.  .l)stop 
k=l 

index=0 
5  iter-iter+1 

if (iter  .gt.  30)  then 
print  1000 
stop 
endif 

16  continue 
zz-zzd 
rz-rzd 
i-1 
ds-20. 

omega-v*l . 69/rpl 
m-1 

arg-1. - .006875*zz/1000. 
if(arg  .gt.  0.)then 

rhoz-0 . 002378*arg**4 .256 
else 
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print  *,  'does  not  converge' 
stop 
end  if 

omegsq=omega* omega 
s=0 . 
thz=0 . 

qs“.  5*rhoz*rz*rz*omegsq*abase 
a=0 . 5*amg*x/qs+acl*x* . 04 
b= ( acl+cddrog) *x - ddrog*acm 
c=-0 . 04*acl*x-amg*x/qs 
alpha=( -b+sqrt(b*b-4 . *a*c) )* . 5/a 
if (alpha  .ge.  .8)  then 
a2=. 85*qs*xcp 
b2=-amg*x- . 85*pi*qs*xcp 
c2=amg*x*pi* . 5- . 85*qs*xcp 
alpha=( -b2-sqrt(b2*b2-4 .*a2*c2) )* . 5/a2 
m“2 

am=amg/g 

endif 

al-- amg/g*rz*omegsq - qs*acl* . 04 

bl=qs* ( ddrog*acm/ (x*cos ( alpha) ) - acl - cddrog) 

cappal=al/bl 

if(cappal  . le .  .04)  then 
trpz=- amg/g*rz*omegsq 
else 

trpz“-amg/g*rz*omegsq+acl*(cappal- . 04)*qs 
endif 

if (alpha  .le.  .04)  then 
tzpz=amg 

trthpz-=0 . 5*rhoz*abase*cddrog*rz*rz*omegsq 
else 

if(m  .eq.  1)  then 

tzpz=arag-acl*(alpha- . 04)*qs 
trthpz“0 . 5*rhoz*abase*cddrog*rz*rz*omegsq 
else 

if(m  .ne.  2)  stop 

tzpz-amg- . 85*qs*sin(alpha)**2 . *cos (alpha) 
trthpz- . 85*qs*sin(alpha)**3 . 
endif 
endif 

tz=sqrt ( trpz*trpz+trthpz* trthpz+tzpz*tzpz ) 

rpz=trpz/tz 

thpz=tr thpz/ ( tz*rz ) 

zpz=tzpz/tz 

rthpz=thpz*rz 

if (mod (index, 3)  .eq.  0  .or.  k  .eq.  4)then 
print  1 

print  21 , s , rz , zz , thz , tz , rpz , zpz , rthpz 
end  if 
8  s-s+ds 
z-zz 
t“tZ 

r-rz 

rp-rpz 

thp-thpz 
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zp— zpz 
kk-1 

10  za-.5*(z+zz) 
ta-. 5*(t+tz) 
ra-. 5*(r+rz) 
thpa=« .  5*  ( thp+thpz ) 
zpa=. 5*(zp+zpz) 
rpa= . 5*(rp+rpz ) 
arg-1 . - . 006875*za/1000 . 
if(arg  .gt.  0.)then 

rho-. 002378*(1 . - . 0068 75*za/1000. )**4 . 256 
else 

print  *,  'does  not  converge' 
stop 
end  if 

thpasq=thpa*thpa 

rasq“ra*ra 

qd“. 5*rho*d*rasq*omegsq*cd 
strth=sqrt(l . -rasq*thpasq) 

trp-=(ta*ra*thpasq-amug/g*ra*omegsq-qd*ra*rpa*thpa*ctrth)*ds+trpz 
trthp=-( - ta*rpa*thpa+qd*strth**3 .+. 5*rho*d*picf*rasq*omegsq)*ds+ 

$  trthpz 

tzp=-(amug-qd*ra*thpa*zpa*strth)*ds+tzpz 

tl-sqrt(trp*trp+trthp*trthp+tzp*tzp) 

rpl-trp/tl 

rl-rz+(rpl+rpz)* . 5*d3 

thpl»trthp/(tl*rl) 

zpl“tzp/tl 

zl=-zz+.  5*(zpl+zpz)*ds 
rthpl=thpl*rl 
thl=thz+( thpz+thpl)* . 5*ds 
if (abs( (zl-z)/z)  .gt.  .001)  go  to  100 
if (abs( (rl-r)/r)  .gt.  .001)  go  to  100 
if (abs( (tl-t)/t)  .gt.  .001)  go  to  100 
32  if(l  .le.  i  .and.  i  .le.  17)  then 
ds“dsstor(i) 
else 

if(i  .eq.  18)  then 

if (s-al+1 . Oe-4  .eq.  0.0)  then 
ds=dds 
else 

if (s-al+1 . Oe-4  .It.  0.0)  then 
go  to  50 
else 

go  to  190 
end  if 
endif 
else 

print  i  is  not  between  1  and  18' 

stop 
endif 
endif 
i-i+1 
50  zz-zl 
rz-rl 
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tz=tl 

rpz=rpl 

thpz=”thpl 

zpz=zpl 

thz=thl 

trpz=trp 

trthpz=trthp 

tzpz^tzp 

if(al-s  .It.  1000 . )ds=al-s 

if(k  .eq.  4)print  21 , s , rl , zl , thl , tl , rpl , zpl , rthpl 
go  to  8 

100  if(kk  .p'.  20)  then 

•  print  3 

go  to  32 

endif 
z=zl 
r=rl 
t=tl 
rp=rpl 
thp=thpl 
zp-zpl 
rthp-rthpl 
kk=kk+l 
go  to  10 
190  continue 

if (mod(index, 3)  .eq.  0  .or.  k  .eq.  4) then 
print  21 , s , rl , zl , thl , tl , rpl , zpl , rthpl 
end  if 

if(iter  .gt.  itrset)go  to  18 
go  to(200,203,18,1199) ,k 
200  if (abs ( (zl-zpl)/zpl) - . 010)202 , 210 , 210 

202  go  to(204,240)  ,k 

203  if (abs (zl- zpl)- 2000. ) 204, 204, 210 

204  if(abs((rl-rpl)/rpl)- .003)206,220,220 
206  go  to(240,200)  ,k 

210  k-1 

zzd-zzd+ . 6*(zpl-zl) 
go  to  5 

220  if(k  .eq.  1)  then 

drzd=. 16*(rpl-rl) 
else 

if(k  .eq.  2)  then 

if (abs ( (rl-rstore)/rpl)  . le .  .002)  drzd=. 2*(rpl-rl) 
.  drzd'=(rpl-  rl)*drzd*.  6/(rl-rstore) 

if(abs(drzd)  .gt.  150.)  then 
drzd=-150 .  *drzd/abs  (drzd) 

•  else 

if (abs (drzd)  .le.  5.)  drzd-5 .*drzd/ abs (drzd) 
endif 
else 

print  k  is  not  1  or  2' 
stop 
endif 
endif 

rzd-rzd+drzd 
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k=2 

rstore-rl 
go  to  5 
240  k-3 

18  continue 

if(index  .gt.  30)go  to  999 
if(istop  .eq.  500)go  to  1199 
if (mod( index . 3)  .eq.  0)then 

if (abs( (zl-zpl)/zpl)  .It.  .001  .and.  abs( (rl-rpl)/rpl)  .It. 

$  .001)go  to  17 

fl=zl-zpl 
f2=rl-rpl 
print  *, ' fl=' , fl 
print  'f2-' , f 2 
zlsav=zl 
rlsav=rl 
zzdsav=z2d 
rzdsav=r2d 
zzd“zzd+delz 
index=index+l 
go  to  16 
else 

if (mod(index, 3)  .eq.  l)then 
dzldzd=(zl-zlsav) /delz 
drld2d= ( rl - rlsav) /delz 
zzd=2zdsav 
rzd=rzdsav+delr 
index“index+l 
go  to  16 
else 

dzldrd=(zl-zlsav) /delr 
drldrd=(rl-rlsav) /delr 
den=dzldzd'*dr  Idrd-  drldzd*dzldrd 
delzd=(f2*dzldrd-fl*drldrd)/den 
delrd= ( f l*dr Idzd- f 2*dzldzd) /den 
zzd=z7:dsav+delzd 
rzd=rzdsav+delrd 
index“index+l 
go  to  16 
end  if 
end  if 
17  k-4 

print  1301 

print  1101,v,rpl,zpl,rzd,zzd 
print  1102, amg, abase, cddrog 
print  1103 , amug, d, al ,dds ,picf ,cd 
print  1203 , alpha, cappal 
istop-500 
go  to  16 
1199  stop 

999  print  *,  'does  not  converge' 
stop 
c 

1  format (Ih  , 5x, Ihs , llx, Ihr , llx, Ihz , llx, 2hth, lOx, Iht , lOx, 2hrp , lOx, 
$  2hzp , 9x,4hrthp) 
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2  format(8el2 .4) 

21  format(8el2 .4) 

3  format(13h  not  converge) 

1000  format ( Ihl , 14hrlnot  converge) 

1101  format (Ih  , 8hv(knot)  =  ,el2 .4, 2x,4hrpl=, el2 .4 , 2x ,4hzpl= , el  2 . 4 , 2x , 

$  4hrzd- , el2 . 4 , 2x,4hzzd“, el2 .4) 

1102  format (Ih  ,4hamg“ ,el2 . 4 ,2x, 6habase=*, el2 .4 , 2x , 7bcddrog= ,  el2 . 4) 

1103  format (Ih  , 5hamug= , el2 .4 , 2x, 2hd=, el2 .4 , 2x, 3hal=, el2 .4 , 2x, 4hdds= , 

$  el2 . 4 , 2x, 5hpicf= , el2 . 4 , 2x, 3hcd=, el2 .4) 

1203  format(lh  , 3x, Shalpha-, el2 . 4 , 5x, 7hcappal= , el2 . 4) 

1301  format (Ihl , 5x, 3hrpl , 9x, Shrdrog, 7x, 3hzpl , 9x, 5hzdrog, 7x, 3hsep , 9x , 

$  Iht , llx , 4hrppl , 8x,4hzppl , 8x, 4hthpl) 

end 
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